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ABSTRACT 


The electric far field radiation pattern is determined for a 
uniformly illuminated, linearly polarized circular aperture 
transmitting through a dielectric radome located in the near 
field of the aperture. A modified electric field integral 
equation is solved using the method of moments procedure and 
the thin shell approximation for dielectrics. The resulting 
solution was computer coded for ogive and spherical radome 
shapes. The program is designed in a modular fashion to 
accommodate the addition of different antenna types, 
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I. INTRODUCTION 

The term radome was coined during World War II as a 
composite of two words, radar and dome. The term originally 
described the class of dome-shaped structures designed to 
house and protect radar antennas on airborne platforms. 
Current usage of the term has evolved to encompass any 
structure that protects or houses an antenna [Ref. 1]. 

Radomes must satisfy a diverse range of engineering 
requirements. These specifications can be broadly divided into 
two categories: electrical and aero-mechanical. Due to aero- 
mechanical considerations (which are not addressed in this 
thesis) an ogive-shaped radome is commonly found on airborne 
platforms. Chapter II discusses the geometry for an ogive- 
shaped body of revolution. Aircraft radomes, especially those 
used at supersonic speeds are subject to mechanical stress and 
aerodynamic heating so severe that the electrical 
characteristics of the radome may be a secondary 
consideration. 

The goal of this thesis is to develop a computer program 
that provides an accurate model for the transmission 
characteristics of a radome, given the geometry and electrical 
characteristics of its material composition. This program is 
called LDBORMM.F, and is described in detail in subsequent 


chapters. 


A. ELECTRICAL CHARACTERISTICS OF RADOMES 

Electromagnetic transmission characteristics are a primary 
electrical design problem for radomes. In this thesis, the 
radome is modeled as a thin isotropic, homogeneous dielectric 
shell. The reflection coefficient (DIT) of the radome is given 


by 


IT = "TX 5 (1.1) 


where R, is the surface resistivity of the radome and n is the 


impedance of free space (1207 ohms). The surface resistance, 


which can be complex in general, is given by 


R, = = (1.2) 
where E, ls the tangential electric field at the surface of 
the radome and J, is the surface current [Ref 2: p. 433]. 
From Equation (1.1), it follows that the reflectivity (T) 
of the radome increases as the surface resistance (R.) 
decreases. In the limit, where R,=0 (perfect conductor), T=-1 
as expected. When an antenna radiates through a radome, a 
small portion of the main beam is reflected off of the walls 
and forms a lobe. This lobe is commonly called the reflection 


lobe or image lobe. The image lobe is usually broader than the 





The primary concern regarding the image lobe is that large 
discrete reflectors in the detection field can become visible 
to the radar through this side lobe, and be interpreted as 


moving targets by the radar [Ref. 3: pp. 301-311]. The reason 


image lobe Y 


scanned 
source 


reflected energy can 
cause significant increase 


` . 
~ main lobe 


N 
in side lobe levels. w A. 





Figure 1.1. Radome Reflection Lobe Geometry. 


for false detection of motion from stationary reflectors is 
the apparent doppler of the return of the image lobe since it 
15 at a different angle and thus a different frequency than 
the main beam lobe. 

The ideal radome is perfectly transparent to the 


electromagnetic energy that must pass through it. The 


The ideal radome is perfectly transparent to the 
electromagnetic energv that must pass through it. The 
materials available and suited for radome construction in high 
velocity airborne systems are not transparent. Impedance 
mismatches (radome/air interface) and radome geometry alter 
the transmission characteristics of the system. The presence 
of a radome can affect the gain, beamwidth, sidelobe level and 
the direction of the boresight, as well as change the VSWR and 


the antenna noise temperature (Ref. 4: p. 265). 


B. METHODS OF ANALYSIS 

Most analytic models for radomes are based on microwave 
optics (ray tracing techniques). These models assume a radome 
curvature large enough so that the surface can be considered 
locally flat. In addition, multiple reflections and surface 
waves are usually ignored. These approximations are not valid 
for systems where the radome is in the near field of the 
antenna or the curvature of the radome causes rapid variation 
in the direction of the tangential vector along the surface of 
the  radome. Consequently, ray tracing techniques are 
inaccurate for the analysis of high performance antenna and 
radome configurations. Figure 1.2 illustrates the inaccuracies 
inherent to ray tracing techniques. Other analytical methods 
are available, but have various shortcomings and limitations. 
Table 1.1 gives a summary of various analytic methods for the 


solution of electromagnetic scattering problems. 


The electric field integral equation (EFIE) and the 


method of moments is employed in this thesis. The objective is 


to obtain the solution for the unknown current density which 
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Figure 1.2. Ray Tracing Approximations. 





occurs in the integral equation. The integral equation is then 
Solved using a numerical procedure called the method of 
moments (MM) [Ref. 5: p. 671]. When properly implemented, a 
method of moments solution is considered rigorous; all of the 
scattering mechanisms are included. This thesis adapts a code 
developed by J.R. Mautz and R.F. Harrington for scattering by 
bodies of revolution. 

The proximity of the antenna to the radome requires a near 
field solution of the integro-differential electric field 
equations. Chapter III develops the analytical solution for 
a radome in the near field of a circular antenna and the 
numerical method for its calculation. 

Chapter IV provides details of the program LDBORMM.F. This 
program solves for the far field radiation pattern of a 
circular aperture transmitting through an arbitrary-shaped 
dielectric body of revolution. The EFIE for electromagnetic 
scattering is solved numerically for the system depicted in 
Figure 1.3. 

Chapter V discusses the testing and evaluation of the 
main program. The effects of radome geometry and dielectric 
properties on the far field radiation pattern are examined. 

The appendices provide copies of the source codes and a 
description of variables in the arguments of the subroutines 


and functions contained in the main program. 


radome 


issile body (body of revolution) 


(conductor) 


— 
boresight 


axis of 
symmetry 


linearly polarized | main beam N 
antenna eleciromicolly 
scanned 





Figure 1.3 Physical Arrangement for LDBORMM.F. 


II. OGIVE GEOMETRY 


A. BACKGROUND 

Forward-looking radar systems in missiles and aircraft are 
by necessity located in the nose of the platform. Aerodynamic 
considerations dictate an ogive-shaped radome for most 
applications. The basic problem is illustrated in Figure 1.3. 
In this thesis, the subroutine OGIVE provides the physical 
dimensions of the surface, where the integro-differential 
field equation solutions are satisfied for a loaded (i.e. 
dielectric) body of revolution (BOR). The subroutine 
TESTSPHERE generates a spherical BOR to compare numerical 
solutions of the main program and function subprograms with 
known analytical solutions. TESTSPHERE is used solely for 
validating the subroutines that calculate the circular 
aperture radiation field. It is not used in the calculation of 
the radome scattered fields. In all configurations the antenna 
is centered on the origin of the global coordinate system 


which is the cylindrical coordinate system of the BOR. 


B. DEFINITION OF OGIVE 

The ogive is a figure of revolution obtained by rotating 
an arc of a circle about an axis in the plane of the arc as 
depicted in Figure 2.1. The curve of the ogive is rotated 


around the z axis to produce a surface of revolution. 


The profile of the surface is generated in the r-z plane 


by the equation 


r(z)s/R*-z'!«b-R (2.1) 


R = radius of parent circle 


OGIVE : body of revolution (bor), r(z) is rotated 
around Z axis The radius (R) of the parent 
circle originates in the 2 * O plane. 





Figure 2.1 Ogive as Defined by Equation (2.1). 


where R is the radius of the parent circle and b is the radius 


of the base. 


10 


C. MODIFIED OGIVE FOR LDBORMM.F 

In the program LDBORMM.F the antenna is located at the 
origin of the coordinate system. In a real system, however, 
the antenna may be displaced on the z axis, from the origin of 
the defining Equation (2.1). The subroutine OGIVE prompts the 
user to enter the value z', which is the displacement on the 
z axis of the antenna location and the position on the z axis 
where the radius (R) of the parent circle of the ogive 


commences rotation. Figure 2.2 depicts the modified profile of 


an ogive as calculated by the subroutine OGIVE. 





Figure 2.2 BOR Generated by Equation (2.2). 


The modified definition for the radius as a function of 


position on the z axis is given by, 


11 


r (2) | j zsz’ (2.2) 
V/R?-(z-z)?*«b-R zz’ 
where z’ is the displacement on the z axis. Using this 
modified definition the source was located at the origin of 
the cylindrical coordinate system used to generate ogive(i.e., 
the source is located in plane z-0). To reiterate, the z' of 
Equation (2.2) and Figure 2.2 denotes the distance from the 


origin where the radius of the parent circle is defined and is 


not related to the source coordinates: 


D. SUBROUTINE OGIVE 

The subroutine OGIVE returns values for the global 
variables NP, ZH and RH. NP is the number of surface 
generating points on the ogive in the r-z plane. In the 
subroutine the range of NP is limited to 3 < NP < 400. The 
program prompts the user for the following inputs: radius of 
the parent circle (R), z', and the base radius (b). 

A standard rule of thumb for the convergence of the series 
approximation for the surface current is that the spacing of 
the surface generating points be approximately one-tenth 
wavelength or less. The following algorithm computes the 
number of points required for a segment size of approximately 


one-tenth of a wavelength on the BOR 


12 


2 -y2*R*b-b? 
9 = ARCSIN (3) 
(9533) 
AL "3 pa Z: 
NP = 2*INTEGER(5*AL)+1. 


NP is the global variable for the number of points in the r-z 
plane of the ogive [reference Figure 2.2]. For a circular 
aperture, the field incident on the radome tends to vary more 
rapidly as the angle off boresight increases, thus more 
generating points are needed. By slightly modifying the above 
algorithm, a surface with any arbitrary segment size can be 
generated. 

Equation (2.3) gives the algorithm employed by the 
subroutine OGIVE to calculate the number of points (NP) on the 
ogive surface in the r-z plane. Data generated uniformly along 
the radome surface will be nonuniform in z, r, and 0. The 
subroutine TESTSPHERE generates data points with a uniform 
spacing in angle 0, whereas OGIVE points are uniform along S. 
Thus a comparison of the two outputs will not yield coincident 
data points. Appendix A contains the source code for 
subroutines OGIVE and TESTSPHERE. Appendix E provides 
mathematical details relevant to the distribution of data 


points in angle 0 for the subroutines OGIVE and TESTSPHERE. 


15 


If the number of points is greater than 400, array limits 
are exceeded and the main program will write a warning to the 
Screen and stop the program. For bodies of revolution where 
the arc length is greater than 40 A the dimensions of the 
following global variables must be changed accordingly 
throughout the program: RH, ZH min R, B, CG ZIO and ZC MEHS 
and RH, are the product of the local coordinates z,, r, 
(computed by the BOR subroutine) and the wave number k=27/A.. 
All spatial dimensions throughout the main progran, 
Subroutines and functions are in wavelengths (A). 

For surface shapes with definitive geometries other than 
an ogive or sphere the user must provide an appropriate 
subroutine. The subroutine must return the global variables 
NP, ZH, RH, b, R and Z’. These variables are required to 
define the generating points of the BOR profile (see Eqn.2.2). 
The call for the subroutine must be inserted in the logical 
block of code in the main program which establishes the BOR 
geometry. The subroutine is appended to the main program. 
Chapter V addresses modification of the main program for 
various bodies of revolution. The program LDBORMM.F saves the 
inputs from the BOR subroutine and writes them to the 


external sequential file OUTLDBOR. 


14 








III. METHOD OF MOMENTS 


A. OVERVIEW 

In order to calculate the scattered field pattern of a 
dielectric radome the electric current (J,) must be integrated 
over the surface of the body of revolution. The magnetic 


vector potential (A) and scalar electric potential (9) are 


defined in terms of the surface current 
E*(R,0,0) * -joA (Jj) - VO(J,) YZ) (3.1) 


where E'(R,0,9) is the scattered electric field in spherical 


coordinates, and L(*) is an operator introduced for notational 


convenience. The potentials are 





p-jkR 
ກບ 2. Aug 9 (3.2) 
5 
bei 
dy fo (3.3) 
E 
a E (3.4) 
jo ວ “+ " 


R is the magnitude of the difference r and r',the position 


vectors of the field and source points respectively. The 
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operator v,* is the surface divergence on S (Ref. 6: pp. 15- 


16]. Figure 3.1 illustrates the quantities. 


P(r,8,0) 


To calculate the electric field , due to a source 
with an electric current density J 
Equations (3.1) through (3.4) must be calculated. 





Figure 3.1. Scattering due to Surface Currents. 


The objective is to solve for the current on the surface 
of the radome. Following the method of moments (MM), the 


electric current (J,) on the surface (S) is represented by 
- t 7t $ 
Js Vi; (Feb Jai * Ie Jai) | (3.5) 


In Equation (3.5) Jt and J* are the  tangential and 
azimuthal components of the current on the surface of the 


radome. The expansion functions for J, are chosen to be 


16 





triangle function 


t y a 


pulse function 


Figure 3.2 Basis Functions Pulse (P) and Triangle (T) of 
Equatlons (5.6) and (3.7) 
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ແ = ejne ດ. , 2... 52 (3.6) 
D 
Pt : 
ປ = d E eine J 71,27 eee 1 (3.7) 
; 
where the azimuthal mode number is n=0,t1,t2....t%. An exact 


solution requires an infinite number of terms in the sum over 
the index n. In practice, when the source (antenna) is on the 
axis of symmetry, the sum converges rapidly and only a few 
azimuthal modes are required. 

The unit vectors t and à are in the tangential and 
azimuthal directions. The func Oneni S ດຍ (t ME 
triangle and pulse functions as shown in Figure 3.2. The 
abscissa t is the arc length along the generating curve of the 
BOR. It is assumed that the generating curve consists of NP-1 
straight line segments where NP is an integer. The expansion 
functions of Equations (3.6) and (3.7) are especially 
appropriate if the BOR is an infinitely thin perfectly 
conducting surface with edges. This is true because the t 
directed electric current approaches zero at an edge, whereas 
the È directed electric current may grow large (Ref. 7: p. 6]. 
The MM solution of the EFIE reduces the integral equation to 


a matrix equation for the unknown coefficients I* and I? . This 
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is accomplished by the MM testing procedure. The test 
functions (VW. Wasa) are chosen to be the complex conjugates of 
Equations (3.6) and (3.7). Equation (3.5) is substituted 
Mto Eguations (3:2) to (3.4) for the «solution of Equation 
(3.1). The dot product of Equation (3.1) is then taken with 
the test functions W,'^. These dot products are then 
integrated over the surface S. The resulting matrix equation 
15 

7 ເເ zt$ -1 
z$t ze 


Tet 
I$ 





E (3.8) 
ve 














where Z is a square matrix called the impedance matrix, and V 
is a column vector called the excitation vector [Ref. 7: pp. 
6-31]. 

The impedance elements are identical to those described by 
Mautz and Harrington. Detailed equations for them appear in 
reference [Ref. 7] and will not be repeated here. However, the 
excitation elements will depend on the characteristics of the 
antenna and the location of the radome. The derivation of the 
excitation elements for a BOR in the near field of a source 
has not been presented elsewhere, and therefore will be 
derived in detail in section (3.B.). 

In LDBORMM.F, the subroutine ZMAT calculates the impedance 
matrix and the subroutine GENEX calculates the excitation 
vector (V). The subroutines DECOMP and SOLVE solve the system 


of equations resulting in the electric current coefficients 
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(I"?). The current coefficients are substituted in Equation 
(3.2), which determines the surface current (J,) [Ref. 7: 
pp. 41-647. 

With the surface current (J,) calculated, the scattered 
electric field (E') is determined by Equation (3.1). The total 


electric field, E(R,0,9), is then given by 


E(R,0,0)- E5 + E? (3.9) 


where E? is the electric field of the circular aperture. 


B. NEAR FIELD OF A CIRCULAR APERTURE 

In order to maintain a clear picture of the analysis that 
follows, it is advantageous to keep in mind the following: 
1. The electric current distribution over the circular 
aperture, which represents an antenna, is defined in Cartesian 
coordinates. 
2. The electric field due to the aperture and the scattered 
electric field due to the radome are computed in spherical 
coordinates (R,9,¢). 
3. Since the ogive is a body of revolution it is generated in 
a cylindrical coordinate system (r,$,z). 
4. Primed quantities denote source coordinates, unprimed 
quantities denote coordinates at the point of observation 
(with the exception of z'/, which was discussed earlier). 

To calculate the excitation vector (V) of Equation (3.8), 


the electric field of the antenna must be known. Since the 
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antenna is typically operating in the near field of the 
radome, far field approximations are not valid and a more 
accurate method of calculation must be employed. Since many of 
the radar antennas used in airborne applications are circular 
paraboloids or slotted waveguide plates, the radiation source 
is modeled as a circular aperture with a uniform current 


distribution polarized in the x direction 


+. Ea (3.10) 
For simplicity the electric current distribution over the 


aperture is assumed to have constant amplitude and a linear 


phase 


J, = Tee EE cos, f (3.11) 


Jo is the amplitude of the electric current, while 0, and 4, 
are scan angles in the spherical coordinate system. The 
coordinates of a point in the plane of the aperture are given 
(X y?! , 0). 

To compute the Cartesian components of the electric field, 
J, is used in Equations (6-108a) to (6-108c) of Reference 
[Ref. 5: p. 283]. A straightforward transformation 61 
rectangular to polar coordinates (see Appendix E) yields the 
following expressions for the Cartesian components of the 


scattered electric field 


2L 








= “IN, x1 A 2 / / / 
E e US [G, *(r cosġ-r' cosġ) G.|I dz AQ (3.12) 
E, = INI ffe” (r cosġ-r' cosġ) 
7 ARK (3413) 


(r sinġ-r'sinġjG,r' drieġ' 


a LE T e'Hr cosġ-r'cosġ)z G, r' dr'dġ' (3.14) 





where 


x1 = jk(r' cosġ'sinb,-R) (3.15) 





R- (is-Tirztrarrisin 9) (3.16) 


k'Rt-i-jkR 
R? 


G1 > (3.17) 


3+3 JkR- k2 R2 
R? 


G2 = (3.18) 


Figure 3.3 defines the guantities used in the eguations. It is 
assumed that the antenna is scanned such that ¢.= 0. 
The rectangular components of the electric field are 


converted to spherical components using the transformation 
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Es = E„cos8cos6 + E,cosBsinė - E,sin8 (3.19) 
Ez = E,Ssin0cos$*E,sinOsin$*E,cosO0 ( 3920) 
E, 8 -E,SinġtE,COBQ . (3.21) 


The function subprograms 
CIRCTHETA , CIRCRHO and CIRCPHI calculate the spherical 
electric field components at the point (P) 


The point (P) is defined in cylindrical coordinates, 
r(2), O and Z 


The primed coordinates are in the plane Z = 0 


The circular aperture Is centered at the origin 





Figure 3.3 Geometry for Equations (3.12) to (3.18). 


The function subprograms CIRCTHETA, CIRCRHO, and CIRCPHI 


calculate the spherical electric field components Eg , E and 


R / 


E,, respectively. 


r4 
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1. FUNCTION SUBPROGRAMS, CIRCTHETA, CIRCRHO, AND CIRCPHI 

The function subprograms CIRCTHETA,  CIRCRHO and 
CIRCPHI are called in the subroutine GENEX. In the main 
program, the vector field due to the source (E?)is added to 
the scattered field vector (E°) of the body of revolution to 
obtain the total field (Ep), per Equation (3-9). Sinse LE 
far field approximation is valid in this case, the closed form 
solution for a circular aperture derived in Section III.B.1.b 
can be used. Subroutine GENEX requires the components of E, in 
the near field. The function subprograms calculate the 
spherical electric field components due to a uniformly 
illuminated circular aperture at an arbitrary point in space 
specified by the spherical coordinates R, 0, $. The source can 
be scanned off boresight by entering non-zero scan angles when 
the program is executed. 

In the subroutine GENEX, these functions are called 
for each point of integration in the tangential and azimuthal 
components of the excitation vector (V) which occur in 
Equation (3.8), and are described in the following section. 

Appendix D contains a table of the variables in the 
argument list of the function subprograms and the source code 
for CIRCSUB.F. The program CIRCSUB.F was used to validate 
the algorithms of  CIRCTHETA, CIRCRHO and CIRCPHI. The 


programs TESTCIRC.M and SCANPLT.M (Appendix D) were used to 
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generate the analytic solution and plot the comparative 
results. 
a. Numerical Method 

To determine the electric field of the antenna and 
the excitation vector requires evaluating several integrals. 
In general these integrals can not be reduced to a closed form 
expression for the geometries under consideration. The 
integrals must therefore be integrated numerically. The chosen 
method for evaluating all integrals in this program is 
Gaussian quadrature. In Gaussian quadrature algorithms, the 
integral is over a normalized interval, such that -1< x < 1. 


N-point Gaussian quadrature approximates an integral by 


IR 


l, fix) Ax Vi n (Aa) (3.22) 


where 
NO ODD (3.23) 


and «9, are the coefficients and x, are the zeros of the 
Legendre polynomial of order m (m-0,41,42...0) [Ref. 8]. 

The number of terms in the sum is alwavs an even 
number. The program GAUS.F (see Appendix C.) will generate an 
external sequential file for any range of N where O < N £ 500. 
The first record in the file is an even number (N). The 
remainder of the file contains N rows of X, and è, . The gaus 


files are stored in a subdirectory named "gaus" of the 
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directory containing the main program LDBORMM.F. The path to 
the data files is "(MAIN DIRECTORY)/gaus/gaus###, where ### is 
N. To create the data file (gaus###), execute GAUS.F in the 
subdirectory "gaus." The screen prompts the user for the data 
filename and an even number of points. 

For integrals that are over a range other than -1 
to +1, a change of variables is necessary before Gaussian 
quadrature is applied. Integrations over the antenna aperture 
(from zero radius to the antenna radius r,) require this 


change of variables is given by [Ref 8: pp. 523-527] 


Di T 
E - (+ l (3.24) 


A similar technique is used to evaluate the ġ 
integrals. The interval of integration for ġ is svmmetric and 


does not require a change of variables since, -T< 0 < T 
O ERX tR . (3.25) 


Equations (3.12) through (3.14) are integrated by this method. 
In the function subprograms, CIRCTHETA, CIRCRHO and CIRCPHI, 
the variables of integration over the aperture are PHIPRIME 
(ġ') and RP (r') and the limits are from 0 to 27 and O to r,, 


respectively. 
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b. Test Cases for Function Subprograms 
In order to verify the algorithms and coding for 
CIRCTHETA, CIRCRHO and CIRCPHI the field magnitude E, was 


computed and compared with the analytical solution given by 


DEAL c I gore 
ະແ ວະ . 
R | 


siné-sin0,) 


( 
| 3.26) 
sin0-sin0,) — 


where E;=nJ, is the electric field amplitude, J, is the first 
order Bessel function, 0, is the scan angle and (R) is the 
distance from the center of the aperture to the far field 
point (P) of Figure 3.3. It is assumed that R is much greater 
than r' in the above formula. Therefore any vector from a 
source point to the observation point is approximately 
parallel to the vector R [Ref. 9: pp. 478-487]. 

The far field distance for which Equation (3.26) is 
valid depends on radius (ນ) of the aperture. Antenna 
engineers commonly consider the far field to exist at 


distances greater than 


R = — (3.27) 





However in optics the far zone applies to values of R where 
the Fresnel number is less than one. Using the latter 
criterion the minimum distance for which Equation (3.26) is 


valid is given by (Ref. 10: pp. 95-96] 
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KEA. 


2 (3.28) 
SL 


mg e 
N 


The program TESTCIRC.M in Appendix D is used to 
calculate the analytical solution for E, in accordance with 
Equation (3.26). The results are compared with the calculated 
field computed by CIRCSUB.F. In Figures (3.4), (3.6), and 
(3.8) the analytical solution is plotted with a solid line 
while the values calculated by the test program  CIRCSUB.F 
are plotted as points. 

Figures (3009 (3.6) and (3.8) show the 
convergence of the numerical solution from CIRCSUB.F as a 
function of the number of integration points with the 
analytical solution given by Equation (3.26) for one, two and 
five wavelength radius circular apertures. Figures (3.5), 
(3.7), and (3.9) show plots of the field magnitudes (E,, E;, 
E,) versus the angle in degrees. Figures 3.10 through 3.12 
show plots from CIRCTHETA and Equation (3.26) for various scan 
angles and aperture radii. (Note that the vertical scale in 
Figures 3.4 to 3.12 is in volts/meter, not dB.) 

Table 3.1 summarizes parameters for Figures (3.4) 
through (3.12). The plotted results verify the excellent 
agreement between the algorithm of the function subprograms 


and the analytical solution. It is important to note the 
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number of points of integration in r” (CNRHO) and q” (CNPHI) 
OTAMEELQEtTONS (3.13) through Equation (3.15) required for 
convergence to the analytical solution of Equation (3.26). 
The number of points of integration in the table are the 
minimum number required for convergence. As the radius of the 
aperture increases, the number of iterations required 
increases as the product of CNRHO and CNPHI. For example, 
Test 3 requires 1200 iterations for convergence to the 


analytic solution. In the program  LDBORMM.F the functions 


SUMMARY FOR TEST CASES. 


2 


"2 
3.10 


subprograms are called five times in the subroutine GENEX 
which is in turn called once for each azimuthal mode. Thus the 
number of calls for the function subprograms is the product of 
the number of integration points in t (NT), number of 
integration points in $ (NPHI), the number of azimuthal modes 


(MODES) and the number of generating points on the surface of 
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the BOR. For a large BOR, the computer run time can become 
significant. It is therefore prudent to test the function 
subprograms for the minimum number of points of integration 
required for convergence before executing the main program. 
2. THE SUBROUTINE GENEX 

The subroutine GENEX calculates the excitation vector 
(V) for Equation (3.8). All variables in the argument list are 
inputs with the exception of the excitation vector. Appendix 
B contains a list of variables in the argument of the 
subroutine GENEX. 


The jth element of V* and V^ are given by 


v? - 2 ff mi, Es ds. j=1,2,...NP-2 (3.29) 
nq ຣ 
vi = a W*.- Es ds, j=1,2,...NP-1 (3.30) 


where the indices j and n are the same as Equations (3.6) and 
(3.7), and W^ is the testing function. The test functions are 
chosen by Galerkin’s method and therefore are the complex 
conjugates of the expansion functions J” [Equations (3.6) and 
(SWALI ES is the pth component electric field of the 
antenna and S denotes the surface of the radome. The 
representation p is a generalized for the spherical components 


R, O and ġ . These components are provided by the subroutines 


CIRCTHETA, CIRCRHO SHE EH 
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Figure 3.4 Test 1. 
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Figure 3.5 Test 1, E(R,0,9) Calculated by CIRCSUB.F. 
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Figure 3.6 Test 2. 
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Figure 3.7 Test 2, E(R,0,4) Calculated by CIRCSUB.F. 
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Figure 3.9 Test 3, E(R,0,9) Calculated by CIRCSUB.F. 
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Figure 3.10 Test 4, 8 Scanned to 30 Degrees. 
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Figure 3.11 Test 5, 8 Scanned to 45 Degrees. 
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Figure 3.12 
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The quantities E,, E, and E, completely define E* and are used 
by GENEX to evaluate Equations (3.28) and (3.29) numerically. 
To further reduce equations (3.28) and (3.29) to a form 
suitable for programming requires evaluation of the dot 
products. Referring to Figure 3.13, the tangent direction 


along the surface can be written as 


Eo Psin h ເ: .? (3.31) 


By the orthogonality properties 


S 
. 0 (3.32) 
t-d 


O O O 


Therefore, from (3.28) and (3.29) v'*=v*=v*=0. The remaining 


elements are 


ລ ES ao. | Ce 

Vai he L, l: CARE Talb "[ T: CO, (Ep (de [as (3.33) 
6 -3nD * - 
vt kf j “IT: SE. EE desf, T SERE: dt dà (3.34) 
| E 

VA): jkedaċ p Mp: 3.35 
Al he IN 3 P. dt dė ( ) 
A; = tia x t; (3.36) 
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——— 


ວ..." :; (3537) 


dq EN 

i-e A. (3.38) 
ເນີດ ກະ (3.39) 
CES COS VAMO). (3.40) 


The quantities are defined in Figures 3.2 and 3.13. The 





Figure 3.13 Geometry for Equations (3.30) through 
(3.39) 


elements of the excitation vector (V) are stored in the global 


variable R, where 
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V,t9+V,R is in R(i) i=1 to NP-2 
Vv, is in R(i+MT) i-1 to NP-1 


and MTENP-2 (Ref. 7:57 pe saw 


C. MODIFICATION OF THE IMPEDANCE MATRIX FOR DIELECTRICS 
1. THIN-SHELL APPROXIMATIONS FOR DIELECTRICS 
For a perfect conductor the surface resistance (R.) is zero. 


To satisfy the boundary condition requires 


E, = -E, (3.41) 


for the tangential components of the incident and scattered 
fields. In this case the incident field is that due to the 
antenna. When a dielectric body is present, the incident 
electric field produces a polarization current rather than a 
conduction current. In this case the boundary conditions state 
that the electric and magnetic fields must be continuous. The 


total electric and magnetic fields are 


(3.42) 


En du 
p: 
E 


Maxwell's equations must be satisfied where 


Vx B= -Jou, ä (3.43) 


and by superposition 
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LE — JIĠ ຍ (3.44) 
So, chat 
Vx H= jwe,£+ jwle-e,) E (3.45) 


where the second term on the right side of Equation (3.44) is 


menzero for e+e, 


The polarization current may be defined as 


=o (e cen rE 1. (3.46) 


Maxwells’ second equation becomes 


Ve = RI (3.47) 


0 


Thus E, and H, can be determined from the polarization current 
ups ra earn free space- For dieleċtric bodies with €,>>1 
or very thin dielectric shells with thickness t «« A, the 
tangential component of the polarization current will be much 
greater than the normal component. In these cases the normal 
component can be neglected and a modified form of the EFIE 


results 


Il 


uo "Fo [p cp KS as 


= E° (3.48) 


where the coefficient of J 


J, is called the surface loading or 


surface impedance (R,). (Frequently it is referred to as a 
surface impedance because it can be complex. Here it will be 


called a resistance to "avoid confusion with the surface 
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impedance approximation in which both electric and magnetic 


currents are allowed.) For a lossless material R is 


S 


imaginary; for a resistive material it is real. 


The new matrix equation corresponding to Equation 


(STIS 
[7] = L ale (3.49) 
where 
Za = ff. 9 ເປ) ds (3.50) 
LE s pn ds . (3.51) 


2. EVALUATING THE IMPEDANCE MATRIX 

The impedance matrix (Z) of Equation (3.8) is 
calculated by the subroutines ZMAT, ZLOAD, and ZTOT. If for 
some reason the impedance of the radome is zero (a perfect 
conductor), only the original impedance matrix from ZMAT is 
needed. For other materials the load impedance is calculated 
by ZLOAD and added to the ZMAT result by the subroutine ZTOT. 

The thin shell approximation assumes the thickness(t) 
is much less than the wavelength (A) and the skin depth (6) of 
the shell material. This is a valid assumption for most 
practical radome designs and materials. The polarization 
current is primarily tangential to the surface of the 
dielectric and the normal components are negligible. This 


approximation is applicable to lossv as well as loss free 
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dielectrics. In the lossy case the load matrix (Z,) has a real 
and imaginary component, while the loss free case is purely 
reactive. 

At this point it is necessary to relate the surface 
impedance to the electrical properties of the material such as 
the dielectric constant and loss tangent. The electric field 
due to a thin dielectric shell is given by equation (3.48) 
with R, defined by 


_ 1 (3.52) 
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where Ae-c-e, (Ref. 11: pp. 531,532]. As stated earlier, the 
dielectric constant of the shell may be complex. Letting 


e=e’-je" gives 





al 
R e a 7 3. 53 
"  jo(t£/-j£"-e,) t 
where £'-£,£,. The electric loss tangent is given by 
// 
ca a or (3.54) 
e d i2. 


where 6 is the loss angle and o=we" is the conductivity 
representing all losses in the medium. Therefore the loss 
tangent is a measure of the power loss in the medium. A medium 
is a good conductor if 0>>we, and a good insulator if o<<we 
(Ref. 2: pp. 342,343]. For we, € 1/60A and thickness t-nA, 


substitution of (3.58) into (3.53) yields 
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- D es 
5s ^ n[j(e,-1) *e,tanó] ' dicc 


Typical radome materials have 2 < e, < 10 and 0.0005 < tan 6 


< 0.005 . Inhomogeneous radomes (i.e. sandwiches) can be 


treated by using an equivalent dielectric constant. 
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IV. PROGRAM DESCRIPTION 


A. DATA FLOW 

The program LDBORMM.F has eight subroutines and four 
function subprograms. Subroutines and functions not covered 
here are documented in reference [Ref. 7]. Table 4.1 gives a 
brief description of all subroutines and function subprograms 


in LDBORMM.F. 


TABLE 4.1 SUMMARY OF SUBPROGRAMS 


OGIVE defines the radome dimensions 


GENEX calculates the excitation vector 


iii EL “27 the moment matrix 


ZLOAD subroutine | calculates the submatrix for the 
dielectric BOR 

ZTOT subroutine | calculates the moment matrix for 
LO BOR 


DECOMP inverts IA moment matrix 


SOLVE subroutine | calculates the current 
coefficients 


PLANE calculate E far field 


BLOG ນ 15 fourth order series expansion of 
nog (xr, fOr x< 0 . 1 . 

CIRCTHETA | functions calculate the p, and $ 

CIRCRHO components of the electric field 

CIRCPHI at a point (P). 
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The block diagram of Figure 4.1 shows the structure of the 


calls for subroutines and functions in the main program. 


IL 
subroutines Sphere or ogive 


enter 






parameters 





load external 
gaus*** 
files 






true 


IMP =0 
alse 


ZTOT 


GENEX 


DECOMP 


SOLVE 






ü È i È 4 i 







OUTLDBOR 

m P . BANGM 
. external BXPOLEM 
Output BCPOLE.M 





Figure 4.1 Structure of main program LDBORMM.F. 
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B. PROGRAM PARAMETERS 
l. NUMBER OF POINTS ON BOR 
The number of points (NP) must be at least 3 but less 
than or equal to 400. This restriction is due to the array 
dimensions set in the code and is not a limitation in the 
solution. The subroutine OGIVE generates segments on the 
surface of the BOR at approximately one tenth of one 
wavelength. Thus, a BOR with an arc length greater than 40 2, 
matrix dimensions are exceeded for the global variables RH, 
ZH, Z, R, B, C, ZLO, and ZL. The subroutine TESTSPHERE allows 
the user to select the number of points on the surface of a 
sphere. The points generated by TESTSPHERE are spaced 
uniformly in the angle 0. 
2. DIRECTORY STRUCTURE 
a. Subdirectory GAUS 
A subdirectory named "gaus" must contain the 
external sequential files with the weights and abscissas for 
the Gaussian quadrature algorithm. The data files are 
generated by the program GAUS.F, which is executed in the 
subdirectory "gaus." Appendix C contains the program GAUS.F. 
The status of the files in the main program is "OLD", thereby 
invoking an error condition if the file does not exist. Figure 


4.2 illustrates the directory structure. 
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MAIN DIRECTORY 






LDBORMM.F = = = = = = = — 


TESTM — — — CTEST.MAT 
h j 
GENPLTM — — execute test then genplt in MATLAB 
environment to generate plots. 











BCPOLE.M 
BXPOLE.M 








GAUS 


CIPOTESI 


ANG.M 
ETHET A.M 
EPHI M 


CIRCSUBF — — — 


TESTCIRCM < ERHO.M 
—- 


CIRCMAT 
SCANPLT.M < ລາຕະ ) 
- PLOTS 





Figure 4.2 Directory Structure. 
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b. External Input/Output Files 


The external sequential data files required for the 
main program are shown in Figure 4.2. The programs that 
generate the files are indicated by dashed lines and arrows. 
The dashed line originates at the generating program and the 
arrow terminates on the file generated. 

c. The Program TEST.M 

The programs TEST.M and TESTCIRC.M must be executed 
before the plot routines GENPLT.M and SCANPLT.M. The program 
TEST.M (or TESTCIRC.M) generates the analytical solution given 
by Equation (3.26). The parameters for BOR radius, antenna 
radius and scan angle are entered by the user . These 
parameters must be the same as the inputs to the main program. 
The routine GENPLT.M plots the analytical solution generated 
by TEST.M and the numerical solution calculated by the main 
program for comparison. A comparison of the two results 
illustrates the effect of the radome on the antenna pattern. 

d. The Subdirectory CIRCTEST 

The subdirectory CIRCTEST contains the programs and 
data files required to test for the number of terms necessary 
for the series approximation of the EFIE to converge to the 
analytic solution. The main program requires extensive 
execution time for large BOR or large antenna. In order to 
minimize run time, the minimum number of integration points 


over the antenna is determined by the execution of the 


sd 


programs in the subdirectory CIRCTEST. The plotting routine 


SCANPLT.M plots the results for comparison. 
Table 4.2 provides a summary of the files in this thesis. 


TABLE 4.2 FILE SUMMARY 
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bal ilu: Mc e S "າ ວມີ nto uv 
` . tete . EM ND > 
sti HOD ACD “+, tun EC M s a e UC M aa et tata 
dr t. ; a 2 See າາາ 

ug : "un EINHEIT eter 
ຈ 
s 


| =" ee cs or euer 
main dir Z main program, solves for 
the radiation pattern of 
radome and antenna 
main dia external data file for 
output of LDBORMM. F 
main dir./ data generated by 
LDBORMM.F for plotting 
wlth GENPLT.M 
main dir./ program to generate 
solution for Eqn.(3.26) 





data generated by TEST.M 
plotting routine 


main dir./ generates gaus### files 
gaus/ 

main dir./ data files for Gaussian 
gaus/ quadrature integration 


CIPCSUDTE main dir./ program to test function 
circtest/ subprograms for minimum 
number of integration 
points 


testcirc.m maln dir P; generates solution for 
CIretest Eqn. (3.26) 

c ir es mat main dir./ data generated by 
circtest TESTCIRC.M 

scanplt.m main dir./ plotting routine 

circtest 


C. MAIN PROGRAM MODIFICATIONS 





The main program is designed in a modular fashion in order 
to facilitate modifications for other radome shapes, antenna 


types or dielectric profiles. 
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1. Modifications for Radome Shapes 

The main program contains subprograms to generate 
ogive-shaped or spherical bodies of revolution. When executed 
the main program prompts the user for the type of BOR desired. 
Additional types of BORs may easily be added to the existing 
menu. The following modifications are required: 
(a) Insert code in main program to write selection number and 
type of shape to screen menu. 
(b) Insert code in logical block that calls subroutine 
generating selected shape. 
(c) Append subroutine program to main program. Subroutine must 
return the variables NP, RH, ZH, b, R and 27 . 
These variables are required for subsequent calculations and 
data entry to external files. 

2. Modifications for Antenna Types 

To modify the main program for an antenna type that 15 
not a circular aperture: 
(a) Substitute the alternate function name in the main program 
where ETF and EPF are assigned. These are the far field 0 and 
ó components of the antenna, so a closed form approximation 
can be used if one exists. 
(b) Substitute the alternate functions in the subroutine 
GENEX, to calculate the variables S1, S2, S3, S4 and S5. 


(C) Append the functions to the main program. 
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To be consistent with the model developed earlier, the 
function subprograms must calculate E(R,0,9$) in spherical 
components at points in the near field. 
3. Modification of Dielectric Profile 

The vector ZLO has a number of elements equal to the 
number of segments (NP-1) on the surface of the loaded body of 
revolution. The elements in the array are the complex values 
of the surface resistance (R,) on each segment of the BOR. As 
written, the main program prompts the user to enter the 
complex surface impedance when the main program is executed. 
As presently coded, the entered value is assigned to each 
segment. In order to examine the effects of a nonuniform 
surface impedance: 
(a) Modify the program to read an external sequential file of 
length (NP-1), with the values of the desired impedance 
profile. 
(b) Store the values in the array 210 in the main program 
before calls to ZLOAD and ZTOT. The order of storage in ZLO 
must correspond to the order of the coordinates stored in RH 


ana ZH. 
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— A 


V. SUMMARY OF RESULTS 


A. VALIDATION OF PROGRAM 

In order to validate and test the main program, the 
surface resistance was assigned a very large value. For very 
large values of surface resistance, the scattered electric 
field (E^) approaches zero. The radiation pattern for a BOR in 
the far field with very large surface impedance is therefore 
the far field pattern of the antenna. For a circular aperture, 
the far field radiation pattern is given by Equation (3.26). 

Figures 5.17 through. 5.3 compare the values for the 
radiation pattern calculated bv the main program and the 
analytic solution generated by TEST.M. The calculated values 
are plotted as (.) and the closed form solution to (3.26) is 
a solid line. The ordinate axis units are volts/meter. The 
abscissa is the angle from the axis of symmetry 0 in degrees. 
The plots of Figures 5.1 to 5.3 show the magnitude of the 
electric field. The actual field values are complex and it is 
the complex values that are used to calculate the excitation 
elements. The figures show the parameters NP, NT, NPHI, CNRHO, 
CNPHI, SCAN, and ARAD as defined in Appendix B. The angle 4 


of the receiver is zero for all test cases. The number of 
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Figure 5.1. Convergence in Far Field, Test 1. 


56 


¿€ = SNIDpeJ a3san3a1oade 


‘bed 0 


Soolbop 


001 06 





qp) 1018 SOT 





(“0““66666) = aouepadut xatTduwuoo 


018 ໄວ 
IHdNO 


aTbue ueos 


Co ASTI 0 





8 = S2UP3SsTp 

€ = LN 

0% = IHdN 

06 = S3UTOd jo iequnu 


* SUALANVUVd 


(da 'WWUOAQT '(PTTOS) TEJTJATEUE 


Sooibop 
0S 





4 WWHOGGTI 30 ISAI 






OS 


001 


OST 


00€ 


osz 


00ເ 


E DOES 125300 


Far Field, Test 2. 


in 


Convergence 


Figure 5.2. 


57 


(0' 066666) = dJuepodur XIaTHAWOD 


G = SNTPpe1 21n319de çZ = 32UP3S1Pp 
OC 4 OHUND C = LN 

Ob ü< IHdNO Ob’ = IHdN 

'beg 0 = aTbue ueos 06 - Squrod JO 13qQUNU 
:SUHLANVAVA 


segalbop 


00ໄ 


qp) LOHd 904 ou 






E UU DEB 30 


( Va Wwaoadd '(prTOS) TeoTAATeue 


So9olbop 


001 06 0 


001 
007 
00ເ 
00) 


E 





009 


J WWHOHG'I JO LSAL 


E dtc: 


Test 3. 


/ 


Figure 5.3. Convergence in Far Field 


58 


points of integration for the antenna were determined by the 
program CIRCSUB.F to insure that aperture  integrations 
converged. Figures 3.4, 3.6, and 3.8 show the results. The 
number of points of integration needed over the BOR surface 
(NT and NPHI) becomes very large as the radial extent (hence 
circumference) of the radome increases. In addition to an 
increase in the number of integration points required, the 
number of azimuthal modes also increases for a large BOR. In 
all test cases the mode number is one (n=+1). Test results 
indicate more points of integration and higher modes are 
required for the calculated values to converge to the analytic 
solution as the antenna radius or sphere radius is increased. 
The execution times were in the range of 2 to 4 hours on a Sun 
Sparcstation 2. 

Figures 5.4 through 5.6 show plots of the magnitude of the 
spherical electric field components for the three test cases. 
As expected the scattered field components ETSCAT and EPSCAT 
are 40 db below the source field (ETF) since in all cases R, 
is 10000. EPF is zero since the antenna radiation is 8 
polarized in this plane. Figures 5.4 through 5.6 are plots of 
field magnitude versus angle O0. 

The plots show the convergence of the calculated values as 
measured against the known analytic solution of a circular 
aperture in the far field. These results verify the program 
for the case where the radome is nearly transparent for large 


real surface resistance in accordance with Equation (1.1). 
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Figure 5.4 Electric Field Components Calculated by LDBORMM.F 
for Test 1. 


60 


100 
100 


- 
< 
O o = o 
ຝງ n tal n 
A. 
m 
n > e N = o uv o ບາ 
2 ; o o 
o o o o ° . : 
o o 
l 
e ° 
o ° 
= 4 
= 
< 
O ° a o 
ຝງ wv (a) un 
+ 
tad 
» 
° = o 
= eo ເລ < ດາ o o o o o o o o 
: o u^ o ທ o n 
o o o ° ດາ N N 4 — 


Figure 5.5 Electric Field Components Calculated by LDBORMM.F 
for Test 2. 
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Figure 5.6 Electric Field Components Calculated by LDBORMM.F 
for Test 3. 
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When the radome is transparent the calculated radiation 


pattern becomes the closed form solution of the far field. 


B. EFFECTS OF RADOME ON TRANSMISSION OF ELECTRIC FIELD 

Test cases were conducted for various radome shapes and 
surface impedances. In Figures 5.7 through 5.9 the calculated 
radiation patterns with the radome are compared to the 
patterns of the isolated antenna. A comparison of the two 
curves illustrates the effect of the radome on the antenna's 
performance. The following parameters were constant for the 
three test cases: 
* resistance = 0+j1700 N; (reflection coefficient T=0.2.), 
* number of azimuthal modes = 2. 
* number of integration points;NT=2,NPHI=40,CNRHO=20,CNPHI=40. 
* antenna radius r, = 2A. 
* radius of the BOR base was 3A (open base-no body structure). 
The radome shapes tested were: 
* TEST 1. a sphere of radius, R=3A4. 
* TEST 2. a cone with the half angle = 30°. 
* TEST 3. an ogive with a parent circle radius R = 7.5A,. 
The forward direction is the region 0 < 90? while the rear 
hemisphere is 8 > 90°. The plotted data in Figures 5.8 and 5.9 
have significant backscatter (02907) due to the shape of the 


radome. The electric field magnitude of the sphere (Figure 
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5.7) for 0290 was below a -20 db threshold. This threshold 
was arbitrarily established for purposes of plotting the data. 
The backscatter from the sphere is negligible since the 
radiation from the antenna at each point on the surface is 
near normal incidence. The cone and the ogive shapes reflect 
energy towards other parts of the surface. These multiple 
reflections raise the sidelobes, especially at wide angles 
from the antenna main beam. 

The coupling of the reflected energy to the antenna is 
assumed to be negligible in this model. This assumption is not 
valid for geometries or dielectric materials which cause 
Significant backscatter, because the reflections from the 
radome perturb the antenna excitation. Not only does this 
effect increase the sidelobes, but also causes an increase in 
the input VSWR. For a well designed radome, this mutual 


interaction is a secondary effect. 


C. RECOMMENDATIONS FOR FURTHER STUDY 

The following suggestions are made for further development 
and improvement of the program LDBORMM.F : 
1. Obtain test data for a real system. Run the program for the 
data and compare the results of the computer model to the 
actual measured test data. 
2. Develop test data from an experimental model constructed in 


the ECE Departmental Transient Electromagnetic Scattering Lab. 
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Compare the measured fields with the predictions of the 
computer model. 

3. Perform a parametric design study to access the impact of 
radome shape, source location and dielectric profiles on the 
radiation pattern. 

4. Perform an exhaustive study of the convergence properties 
of the MM solution. 

5. To improve the execution time some of the Gaussian 
integration loops can be replaced with a delta function 
approximation for the integral. 

6. Symmetry between the positive and negative modes can also 
be exploited to reduce run time. 

This thesis has developed flexible, modular code for 
modeling the electric field radiation pattern for a source 
enclosed in a dielectric body of revolution. Further work 
should concentrate on validation and improvements to reduce 
execution time. After these enhancements are incorporated, the 


code can be applied to the radome deSign process. 
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Figure 5.7 Effects of Radome on Far Field Radiation 
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Figure 5.8 Effects of Radome on Far Field Radiation 
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Figure 5.9 Effects of Radome on Far Field Radiation 
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APPENDIX A. SOURCE CODE AND DATA FILE 


SOURCE CODE FOR MAIN PROGRAM 


PROGRAM :  ldbormm.f 
DATE + 23 January 1992 
REVISED : 21 February 1992 
PROGRAMMERS : D. JENN, R. FRANCIS 
>> SPECIALIZED FOR ARBITRARY CIRCULAR APERTURE EXCITATION << 
BE C ITATION IS SPECIFIED INTHE SUBROUTINE GENEX (...-) << 
BASED ON MAUTZ AND HARRINGTON'S COMPUTER PROG FOR BORS. 
all or part of the surface mav be covered with a surface 
impedance. 
SOURCE IS AT BOR COORDINATE SYSTEM ORIGIN. 
imp=0 perfect conductor 
iprint=0 print pattern points to unit 8 
iseg=0 print the generating curve points to unit 8 


BEGIN MAIN PROGRAM*XXXXXXXXXX%X 


CHARACTER*8 GNT,GNPHI,CGR,CGP 

CHARACTER*14 TPTS,PPTS,PHIPTS,RHOPTS 

COMPLEX EP,ET,Z(100000),R(1600),B(1600),C(800),U,UC 
COMPLEXETI,EP1,ET2,EP2,EC,EX,ZLO(400),ZL(2400),ETF,EPF 
COMPLEX EXP1,EXP2,CONJG,CEXP,CMPLX,CT(3000),IMP,JK 
COMPLEX CIRCTHETA,CIRCPHI 

DIMENSION RH(400),ZH(400),XT(4) , AT(4) , IPS(800) 
DIMENSION A(100),X(100),EXP(500) , ANG(500) , ECP(500) 
DIMENSION XP(100),AP(100),XR(50) , AR(50) 

DIMENSION ECV(500) , EXV(500) , PHC(500) , PHX(500) 

REAL ETSCAT (500) , EPSCAT (500) ,ETHF (500) ,EPHF (500) 
INTEGER NT,NPHI,CNRHO,CNPHI,NP,SELECTION 

REAL MODE,BASE,RS,ZP,RHB,ZHB 

DATA PI,START,STOP/3.1415926,0.,90./ 

DATA IPRINT/0/ 

Rad-PI/180. 

ECX-0. 

BK-2.*PI 

U=(0.,1.) 
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Uo dO 
UC-cU UAE 
JK=(070,6. 283265307) 
C Select the subroutine to generate BOR. 
WRITE(6,*) ’MAKE BOR SELECTION’ 


WRITE (16.30% f 
WRITE(6,*)'Enter number from menu to make selection' 
WRITE(6,X)' 2 
WRITE(6,*)"' Selection Number BOR Geometry' 
WRITE(6,*)' 1 OGIVE' 
WRITE(6,*)' 2 SPHERE’ 


READ(5,*) SELECTION 
IF (SELECTION.EQ.1) THEN 
CALL OGIVE(NP, ZH, RH, BASE,RS, ZP) 
ELSEIF (SELECTION.EQ.2) THEN 
CALL TESTSPHERE (NP, ZH, RH, BASE,RS, ZP) 
ENDIF 


CXXXXXXXXXCALL OGIVE or TESTSPHERE. 
IF (NP.GT.399)THEN 
WRITE (6,*) “MAXIMUM NUMBER OF POINTS(NP) IS 399’ 
GOTO 999 
ENDIF 
DT-90.0/FLOAT (NP-1) 
WRITE(6,*)'ENTER THE FILENAMES gaus###’ 
WRITE(6,*)*'....... mau 4“, 
WRITE(6,*)'ENTER THE FILENAME IN T (NT)'” 
READ(5,*)GNT 
WRITE(6,*)'ENTER THE FILENAME IN PHI (NPHI) ’ 
READ(5,*)GNPHI 
WRITE(6,*)'ENTER THE FILENAME IN CNRHO' 
READ(5,*)CGR 
WRITE(6,*)'ENTER THE FILENAME IN CNPHI' 
READ(5,*)CGP 
C OPEN THE FILES FOR THE gaus/gaus### 
TPTS-'/gaus/'//GNT 
PPTS-'gaus/'//GNPHI 
RHOPTS-'gaus/'//CGR 
PHIPTS='gaus/'//CGP 
OPEN(1,FILE=TPTS,STATUS=""" 
OPEN (2,FILE=PPTS,STANM 
OPEN(3,FILE=RHOPTS,! ) 
OPEN (4, FILE=PHIPTS ,STASUe™ VLD: ) 
READ(1,*)NT 
IF(NT.GT.4)THEN 
WRITE(6,*)'MAXIMUM NUMBER OF POINTS(NT) IS 4” 
GOTO 999 
ENDIF 
READ (2, *) NPHI 
IF(NPHI.GT.200)THEN 
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WRITE (6,*) "MAXIMUM NUMBER OF POINTS(NPHI) IS 200“ 
GOTO 999 
ENDIF 
READ (3, *) CNRHO 
IF(CNRHO.GT.50)THEN 
WRITE(6,*) “MAXIMUM NUMBER OF POINTS(CNRHO) IS 50’ 
GOTO 999 
ENDIF 
READ (4, *) CNPHI 
IF (CNPHI.GT.100) THEN 
WRITE (6,*) “MAXIMUM NUMBER OF POINTS(CNPHI) IS 100’ 
GOTO 999 
ENDIF 


C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 


(O O Oy O 


DO 1 M-1,NT 
READ(1,*,END-21)XT(M) , AT(M) 

CONTINUE 

DO 2 M-1,NPHI 
READ(2,*,END-2) X(M) , A(M) 

CONTINUE 

DO 3 M=1, CNRHO 
READ(3,*,END=3)XR(M) ,AR(M) 

CONTINUE 

DO 4 M=1,CNPHI 
READ(4,*,END=4)XP(M) ,AP(M) 

CONTINUE 

GOS (1) 

CLOSE (2) 

CLOSE (3) 

CLOSE (4) 

MP-NP-1 

MT-MP-1 

N=MT+MP 


WRITES, *) “ENTER MODE (FLOAT)' 

READ(5,*) MODE 

WRITE(6,*)’ENTER PHI (observation) IN DEGREES’ 
READ (5,*)P 

PHI=P*RAD 

WRITE(6,*) ‘ENTER THE SCAN ANGLE IN DEGREES’ 
READS, *)SA 

SCAN--SA*RAD 

WRITE(6,*)'ENTER COMPLEX IMPEDANCE' 
READ(5,*) IMP 

WRITE(6,*) ’ENTER THE ANTENNA RADIUS (wavelengths) ’ 
READ (5, *) ARAD 


ENTER THE NUMBER OF AZIMUTHAL MODES 


(ns-MODES,...,0,...,4MODES) 


MODES-MODE 


dl 


NBLOCK=2*MODES+1 
MHI=MODES+1 
OPEN(8,FILE-'outldbor') 
WRITE(8,8000) 2.*BASE,NP,RS,ZP 
8000 FORMAT(//,5X,'*** BOR RADIATION PATTERN FOR A CIRCULAR' 
*' DISC RADIATOR USING GENEX **x'/, 
*//,2X,'BOR DIAMETER (WAVELENGTHS)-',F5.2,/,2X, 
* ‘NUMBER OF GENERATING POINTS (NP)=’,14, 
*//,2X,'SURFACE RADIUS'/,F5.2,2X,'/ZPRIME',F5.2) 
WRITE(8,30) NT,NPHI 
30 FORMAT(/,12X," NT  NPHI',/,10X,15,2X,I5) 
IF(ISEG.EQ.0) WRITE(8,1300) 
1300 FORMAT(/,10X, 'INDEX',8X, 'Z(I)',10X,'RHO(I)',12X,'SURF 
IMPED') 
DO 52 I-1,NP 
IF(ABS(ZH(I)).LT..001) ZH(I) 
IF(ABS(RH(I)).LT..001) RH(I) 
ZHB-ZH(I)/BK 
RHB-RH(I)/BK 
Z10(i)=IMP 


0 
0 . 


C 
C ASSIGN SURFACE IMPEDANCE AT THIS POINT. THE SURFACE 
C IMPEDANCE OF SEGMENT I IS ZCO(T) 


C 
WRITE(8,8004) I,ZHB,RHB,ZLO(I) 

52 CONTINUE 

8004 FORMAT(11X,14,4X,F8.3,8X,F8.3,6X,2F8.2) 

e 


C MODE LOOP TO CALCULATE THE CURRENT COEFFICIENTS. POS AND NEG 
C MODES DONE IN THE SAME ITERATION OF THE LOOP 


Ckkkkkkk**k****7LOAD,ZMAT,GENEX,DECOMP,SOLVE 
IF(CABS(IMP).NE.0) CALL ZLOAD(NP,RH,ZH,ZLO,ZL) 
DO 400 M-1,MHI 
NM-M-1 
CALL ZMAT(NM,NM,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 
IF(CABS(IMP).NE.O) CALL ZTOT(MT,MP,ZL,Z) 
CALL GENEX(NM,NP,NT,NPHI,CNRHO,CNPHI,XT,AT,X,A, 
* XR,AR, XP,AP,SCAN, PHI, ARAD,RH, ZH, B) 
CALL DECOMP(N,IPS,Z) 
CALL SOLVE(N,IPS,Z,B,C) 


Ckkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
C STORE CURRENT COEFFICIENTS IN ONE LONG COLUMN VECTOR 
NTOP1-MODES-NM 
NTOP2=NBLOCK- (NTOP1+1) 
NS2=NTOP1*N 
NS1=NTOP2*N 
C POSITIVE MODE 
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DO 401 L-1,N 
401 CT(L+NS1)=C(L) 
IF(NM.NE.0) THEN 
C NEGATIVE MODE 
NMN=-NM 


CKHEKKKKKKKKKKKKKKKKKKKKZMAT,GENEX, DECOMP,SOLVE. 
CALL ZMAT(NMN,NMN,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 
IF(CABS(IMP).NE.0) CALL ZTOT(MT,MP,ZL,Z) 
CALL GENEX(NMN,NP,NT,NPHI,CNRHO,CNPHI,XT,AT,X,A, 
* XR,AR,XP,AP,SCAN,PHI,ARAD,RH,ZH,B) 
CALL DECOMP(N,IPS,Z) 
CALL SOLVE(N,IPS,Z,B,C) 

C * k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K 


DO 402 L=1,N 
402 CT(L+NS2)=C(L) 
ENDIF 
400 CONTINUE 
IT=NP 
DT=STOP/FLOAT(NP-1.) 
C 
C BEGIN PATTERN LOOP FROM START TO STOP IN INCREMENTS OF DT 
© (ALL IN. DEG) 
ດ 
DO 500 I-1,IT 
THETA=FLOAT(I-1)*DT+START 
THX=THETA*RAD 
PHR=PHRO 
RHB=RH (1) /BK 
ZHB=ZH (1) /BK 
IF(THETA.GT.180.) THEN 
PHR=PHRO+PI 
THX=(360.-THETA) *RAD 
ENDIF 
ET1=(0.,0.) 
EP1=(0.,0.) 
ET2=(0.,0.) 
EP2=(0.,0.) 
DO 300 M=1,MHI 
NM=M-1 
EXP1=CEXP(CMPLX(0.,FLOAT(NM) *PHR) ) 
EXP2-CONJG (EXP1) 


Qc kk kc k kk kk kkkkkkkkXkPLANE 
CALL PLANE(NM,NM,NP,NT,RH,ZH,XT,AT,THX,R) 
NTOP1-MODES-NM 
NTOP2=NBLOCK- (NTOP1+1) 
NS2=NTOP1*N 


73 


NS1-NTOP2*N 
DO 250 L-1,MT 
ET1=ET1+R(L) *CT(L+NS1) *EXP1 
EP1=EP1+R(L+N) *CT(L+NS1) *EXP1 
IF(NM.EQ.0) GO TO 250 
ET1=ET1+R(L) *CT(L+NS2) *EXP2 
EP1=EP1-R(L+N) *CT(L+NS2) *EXP2 

250 CONTINUE 
DO 260 L=1,MP 
ET2=ET2+R(L+MT) *CT(L+NS1+MT) *EXP1 
EP2=EP2+R(L+MT+N) *CT(L+NS1+MT) *EXP1 
IF(NM.EQ.O) GO TO 260 
ET2=ET2-R(L+MT) *CT (L+NS2+MT) *EXP2 
EP2=EP2+R(L+MT+N) *CT(L+NS2+MT) *EXP2 

260 CONTINUE 

300 CONTINUE 

C DISK CONTRIBUTION IN THE FAR FIELD IS ETF,EPF 
RRR=SQRT (RHB**2+ZHB**2) 


Ckkkkkkkkkk****CIRCTHETA,CIRCPHI 
ETF=CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 


* SCAN, PHI, RHB, ZHB) 
EPF=CIRCPHI (CNPHI, XP, AP, CNRHO, XR,AR,ARAD, 
* SCAN, PHI, RHB, ZHB) 


c i deleted the rrr*cexp(jk*rrr) factor from ETF & EPF. 
C * xk k xk k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 


ETHF(I)=CABS(ETF) 
EPHF(I)=CABS(EPF) 
C TOTAL E-THETA AND E-PHI COMPONENTS 

ET=ET1+ET2+ETF 
EP=EP1+EP2+EPF 
ETSCAT(I)=CABS(ET-ETF) 
EPSCAT(I)=CABS(EP-EPF) 
EC=ET 
EX=EP 
ECV(I)=CABS(EC) 
EXV(I)=CABS(EX) 
ECR=REAL(EC) 
ECI=AIMAG (EC) 
EXR=REAL(EX) 
EXI=AIMAG(EX) 
PHC(I)=ATAN2(ECI,ECR+1.e-10)/RAD 
PHX(I)=ATAN2(EXI,EXR+1.e-10)/RAD 
ANG(I)=THETA 
ECX=AMAX1 (ECX,ECV(I) ,EXV(I)) 

500 CONTINUE 
WRITE(6,*) “MAX E VALUE-',ECX 
WRITE(8,103) P,ECX 

103 FORMAT(/,10X,'PHI OF RECEIVER (DEG)-',F8.2,/,10X, 
* 'MAXIMUM FIELD VALUE (V/M)=',E15.5) 


74 


600 


5015 


5016 
9000 


DO 600 I=1,IT 
ECP(I)=(ECV(I)/ECX)**2 
EXP(I)=(EXV(I)/ECX)**2 
ECP(I)=AMAX1(ECP(I),1.E-10) 
EXP(I)=AMAX1(EXP(I),1.E-10) 
ECP(I)=10.*ALOG10(ECP(I)) 
EXP(I)=10.*ALOG10(EXP(I)) 
CONTINUE 

IF(IPRINT.EQ.0) THEN 
WRITE(8,5015) 
FORMAT(//,7X,'ANGLE'/,15X,'CO-POL'/,25X, ' X- POL' ,/ ,7X, 


1'(DEG)/,4X,2(' (VOLTS) ^ , 4X," (DEG) ^ , 3X, ' (DB-REL) / , 4X)) 


DO 9000 L-1,IT 
WRITE(8,5016)ANG(L),ECV(L),PHC(L),ECP(L) , EXV(L) , PHX(L) 


1,EXP(L) 


FORMAT(5X,F6.2,3X,2(F8.4,3X,F7.2,3X,F7.2,3X)) 
CONTINUE 

ENDIF 
OPEN(2,file-'bang.m') 
OPEN(3,file—'bcpole.m') 
OPEN (4,file=’bxpole.m’) 
OPEN (7, file=’etscat.m’) 
OPEN(8,file=’epscat.m’) 
OPEN (9, file=’etf.m’) 
OPEN(10,file=’epf.m’) 
DO 9097 I=1,IT 
WRITE(2,5019) ANG(I) 
WRITE(3,5019) ECV(I) 
WRITE (4,5019) EXV(I) 


C change ecv,ecv to ecp & exp to get normalized db values. 


2097 
5019 


299 


WRITE(7,5019) ETSCAT(I) 
WRITE(8,5019) EPSCAT(I) 
WRITE(9,5019) ETHF(I) 
WRITE(10,5019) EPHF(I) 
FORMAT(F8.3) 

CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

CLOSE(7) 

CLOSE (8) 

CLOSE (9) 

CLOSE (10) 

CONTINUE 

STOP 

END 


CKXXXXXXXXXXX&XXEND OF MAIN PROGRAM. 
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e 


C*k k kk kkkkkkkk**SUBROUTINE ZMAT. 
C REFERENCE: AN IMPROVED E-FIELD SOLUTION FOR CONDUCTING BOR 


() () GO) 0) 


C) Cy C O 


J.R.MAUTZ AND R.F. HARRINGTON 
TECHNICAL REPORT TR-80-1 

ROME AIR DEVELOPMENT CENTER 
CONTRACT NO. F-<-30602-79=C-0011 


SUBROUTINE ZMAT(M1,M2,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 


COMPUTE THE MM IMPEDANCE MATRIX ELEMENTS. 
MAUTZ AND HARRINGTON (NO CHANGES). 


THIS IS FROM 


COMPLEXZ (100000) ,U1,U2,U3,U4,U5,U6,U7,U8,U9,UA,UB,G4A(4) ,G5A 


(2) 


10 


COMPLEXG6A(4),G4B(4),G5B(4),G6B(4) H4A,H5A,H6A,H4B,H5B 


COMPLEX H6B,UC,UD,GA(100) ,GB(100) 


DIMENSION RH(400),ZH(400),X(100),A(100),XT(4),AT(4) 
DIMENSION RS(400) ,ZS(400) ,D(400) ,DR(400) ,DZ (400) 
DIMENSION DM(400) ,C2(100) ,C3(100) ,R2(4) ,Z2(4) 


CT=2. 

CP=.1 

DO 10 I=2,NP 

I2=I-1 
RS(I2)=.5*(RH(I)+RH(I2)) 
ZS(I2)=.5*(ZH(I)+ZH(I2)) 
D1-.5*(RH(I)-RH(I2)) 
D2-.5*(ZH(I)-ZH(I2)) 
D(12)=SORT(D1*D1+D2*D2) 
DR(12)=D1 

DZ (12) =D2 


IF(RS(I2).EQ.0.) RS(I2)-1. 


DM(I2)-2D(I2)/RS(I2) 
CONTINUE 
M3=M2-M1+1 
M4=M1-1 
PI2=1.570796 

DO 11 K=1,NPHI 
PH=PI2*(X(K)+1.) 
C2 (K) =PH*PH 
SN-SIN(.5*PH) 

C3 (K)=4.*SN*SN 
A1-PI2*A(K) 
D4=.5*A1*C3(K) 
D5=A1*COS (PH) 
D6=A1*SIN(PH) 
M5=K 

DO 29 M=1,M3 
PHM- (M44M) *PH 
A2=COS (PHM) 
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DIMENSION C4(100),C5(100),C6(100),27(4),R7(4),28(4) , R8(4) 


23 
TT 


12 


26 


C4 (M5)=D4*A2 
C5(M5)=D5*A2 

C6 (M5) =D6*SIN (PHM) 
M5=M5+NPHI 
CONTINUE 

CONTINUE 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2N=MT*N 

N2=N*N 

U1=(0.,.5) 
U2=(0.,2.) 
JN=-1-N 

DO 15 JQ-1,MP 
KO=2 

IF(JQ.EQ.1) KQ-1 
IF(JQ.EQ.MP) KQ=3 
R1=RS (JQ) 

Z1=ZS (JQ) 
D1-D(JQ) 
D2-DR(JQ) 

D3=DZ (JQ) 
D4=D2/R1 

D5=DM (JO) 
SV=D2/D1 

CV=D3/D1 

T6=CT*D1 
T62=T6+D1 
T62=T62*T62 
R6=CP*R1 
R62=R6*R6 

DO 12 L=1,NT 
R2(L)=R1+D2*XT(L) 
Z2(L)=Z1+D3*XT(L) 
CONTINUE 
U3=D2*U1 

U4=D3*U1 

DO 16 IP=1,MP 
R3=RS (IP) 

Z3=ZS (IP) 
R4=R1-R3 

24-21-23 
FM=R4*SV+Z4*CV 
PHM=ABS (FM) 
PH=ABS (R4*CV-Z4*SV) 
D6=PH 
IF(PHM.LE.D1) GO TO 26 
D6=PHM-D1 

D6=SORT (D6*D6+PH* PH) 
IF(IP.EQ.JQ) GO TO 27 
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40 


35 


34 


37 


38 


KP=1 
IF(T6.GT.D6) KP=2 
IF(R6.GT.D6) KP=3 

60 10 28 

KP=4 

GO TO (41,42,41,42),KP 
DO 40 L-1,NT 

D7=R2 (L) -R3 

D8=Z2 (L)-Z3 
Z7(L)=D7*D7+D8*D8 

R7 (L) =R3*R2 (L) 
Z8(L)=.25*Z7(L) 
R8(L)=.25*R7(L) 
CONTINUE 
Z4=R4*R4+Z4*Z4 
R4=R3*R1 

R5=.5*R3*SV 

DO 33 K=1,NPHI 

A1=C3 (K) 

RR=Z4+R4*A1 

UA=0. 

UB=0. 

IF(RR.LT.T62) GO TO 34 
DO 35 L=1,NT 
R=SQRT(Z7(L)+R7(L)*A1) 
SN=-SIN(R) 

CS=COS (R) 
UC=AT(L)/R*CMPLX(CS,SN) 
UA=UA+UC 
UB=XT (L) *UC+UB 
CONTINUE 

GO TO 36 

DO 37 L-1,NT 
R-SQRT (Z8 (L) *R8(L) *A1) 
SN--SIN(R) 

CS=COS (R) 
UC=AT(L)/R*SN*CMPLX(-SN,CS) 
UA=UA+UC 
UB=XT (L) *UC+UB 
CONTINUE 

A2=FM+RS*A1 
D9=RR-A2*A2 

R=ABS (A2) 

D7=R-D1 

D8=R+D1 

D6=SORT (D8*D8+D9) 
R=SQRT (D7*D7+D9) 
IF(D7.GE.0.) GO TO 38 
A1=ALOG ( (D8+D6) * (-D7+R)/D9)/D1 
GO TO 39 
A1=ALOG((D8+D6)/(D7+R))/D1 


78 


22 


36 


33 


46 


45 


65 


64 


66 
63 


UA=A1+UA 
UB=A2*(4./(D6+R)-A1)/D1+UB 
GA(K)=UA 
GB(K)=UB 
CONTINUE 

K1=0 

DO 45 M=1,M3 
H4A=0. 

H5A-0. 

H6A-0. 

H4B=0. 

H5B-0. 

H6B-0. 

DO 46 K-1,NPHI 
K1=K1+1 

D6=C4 (K1) 

D7=C5 (K1) 

D8=C6 (K1) 

UA=GA (K) 

UB=GB (K) 
H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
H5B=D7 *UB+H5B 
H6B=D8*UB+H6B 
CONTINUE 
G4A(M)=H4A 
G5A(M)=H5A 
G6A(M)=H6A 
G4B(M)=H4B 
G5B(M)=H5B 
G6B(M)-H6B 
CONTINUE 
IF(KP.NE.4) GO TO 47 
A2-D1/(PI2*R1) 
D6-2./D1 

D8-0. 

DO 63 K-1,NPHI 
A1=R4*C2(K) 
R=R4*C3 (K) 
IER LT. T62) GO TO 64 
D7=0. 

DO 65 L=1,NT 
D7=D7+AT(L)/SQRT(Z7(L)+A1) 
CONTINUE 

GO TO 66 
A1=A2/(X(K)+1.) 
D7=D6*ALOG (A1+S0RT(1.+A1*A1)) 
D8=D8+A(K) *D7 
CONTINUE 


T9 
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41 


25 


17 


62 


51 


32 


A1-.5*A2 
A2=1./A1 
D8=-PI2*D8+2./R1*(BLOG(A2)+A2*BLOG(A1)) 
DO 67 M=1,M3 

G5A (M) =D8+G5A (M) 
CONTINUE 

GO TO 47 

DO 25 M=1,M3 

G4A(M)=0. 

G5A(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)-0. 

G6B(M)=0. 

CONTINUE 

DO 13 L=1,NT 

A1=R2(L) 

R4=A1-R3 

74-22(L)-23 
Z4=R4*R4+Z4*Z4 
R4=R3*A1 

DO 17 K=1,NPHI 
R=SQRT(Z4+R4*C3(K)) 
SN--SIN(R) 

CS=COS (R) 
GA(K)=CMPLX(CS,SN)/R 
CONTINUE 

D6-0. 

IF(R62.LE.Z4) GO TO 51 
DO 62 K-1,NPHI 
D6=D6+A(K)/SQRT(Z4+R4*C2(K)) 
CONTINUE 
Z4=3.141593/SORT(Z4/R4) 
D6=-PI2*D6+ALOG (Z4+SQRT(1.+Z4*Z4))/SORT(R4) 
A1=AT (L) 

A2=XT(L)*A1 

K1=0 

DO 30 M=1,M3 

8550 . 

U6-0. 

U7-0. 

DO 32 K-1,NPHI 
UA-GA(K) 

K1=K1+1 

U5=C4 (K1) *UA+U5 

U6=C5 (K1) *UA+U6 

U7=C6 (K1) *UA+U7 
CONTINUE 

U6=D6+U6 
G4A(M)=A1*U5+G4A (M) 
G5A(M)=A1*U6+G5A (M) 
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18 


21 


19 


G6A(M)=A1*U7+G6A(M) 
G4B(M)=A2*U5+G4B(M) 
G5B(M)=A2*U6+G5B(M) 
G6B(M)=A2*U7+G6B(M) 
CONTINUE 

CONTINUE 

A1=DR(IP) 

UA-A1*U3 

UB=DZ (IP) *U4 
A2-D(IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JM-JN 

DO 31 M-1,M3 
FM=M4+M 
A1=FM*DM(IP) 
H5A-G5A (M) 
H5B=G5B(M) 
H4A=G4A(M)+H5A 
H4B=G4B(M)+H5B 
H6A=G6A(M) 
H6B=G6B(M) 
U7=UA*H5A+UB*H4A 
U8=UA*H5B+UB*H4B 
U5=U7-U8 

U6=U7+U8 
U7=-U1*H4A 
U8=D6*H6A 
U9=D6*H6B-A1*H4A 
UC=D7* (H6A+D4 *H6B) 
UD=FM*D5*H4A 
K1=IP+JM 

K2=K1+1 

K3=K1+N 

K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 
K8=K4+N2N 

GO TO (18,20,19),KQ 
Z (K6) =U8+U9 
BBIP.BO.1) GO TO 21 
Z(K3)=Z (K3)+U6-U7 
Z(K7)=Z (K7)+UC-UD 
MRIEPESEHO.MP) GO TO 22 
Z(K4)=U6+U7 
Z(K8)=UC+UD 

CORTO 22 
Z(K5)=Z(K5)+U8-U9 
MUN SEO.1)" GO TO 23 
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29 


20 


24 


Ze 


Su 


13 


Z(K1)=Z(K1)+U5+U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 
Z(K2)=Z(K2)+U5-U7 
Z(K8)=UC+UD 

GO TO 22 

Z(K5)=Z(K5) +U8-U9 
Z(K6) =U8+U9 
IF(IP.EQ.1) GO TO 24 
Z (K1)=Z(K1)+U5+U7 
Z(K3)=Z(K3) +U6-U7 
Z(K7)=Z(K7) +UC-UD 
IF(IP.EQ.MP) GO TO 22 
z (K2)=Z (K2)+U5-U7 
Z(K4)=U6+U7 

Z (K8) =UC+UD 

Z(K8+MT) =U2* (D8* (H52 «H5B) -A1*UD) 
JM=JM+N2 

CONTINUE 

CONTINUE 

IN=JIN+N 

CONTINUE 

RETURN 

END 


CAXKEKEAAAS AX AA SUBROULING SOLVE. 
C REFERENCE : TECHNICAL FEPORT TR-20-1 


G 
C SEE 
E 


SUBROUTINE SOLVE(N,IPS,UL,B,X) 
MAUTZ £ HARRINGTON FOR DETAILS 


COMPLEX UL(100000),B(800),X(800),SUM 
DIMENSION IPS(800) 
NP=N+1 

IP=IPS(1) 
X(1)=B(IP) 

DO 2 I=2,N 
IP-IPS(I) 

IPB-IP 

IM1=I-1 

SUM=(0.,0.) 

DO 1 J=1,IM1 
SUM=SUM+UL(IP) *X(J) 
IP=IP+N 

X(I)=B(IPB) -SUM 
K2=N* (N-1) 

IP=IPS (N)+K2 
X(N)-X(N)/UL(IP) 

DO 4 IBACK-2,N 
I-NP-IBACK 
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K2-K2-N 
IPI=IPS(I)+K2 
IPl-IH1 

SUM-(0.,0.) 

IP-IPI 

DO 3 J=IP1,N 
IP=IP+N 
SUM=SUM+UL(1IP)*X(J) 
X(I)-(X(I)-SUM)/UL(IPI) 
RETURN 

END 


C******SUBROUTINE DECOMP 


C 


E 


EEEERENCE: TECHNICAL REPORT IR-80-1 


SUBROUTINE DECOMP(N,IPS,UL) 


C SEE MAUTZ & HARRINGTON FOR DETAILS 


E 


Np 


10 


11 


14 


15 


COMPLEX UL(100000),PIVOT,EM 
DIMENSION SCL(800),IPS(800) 
DO 5 I-1,N 

IPS(I)=I 

RN=0. 

J1=I 

Do 2 J=1,N 

ULM=ABS (REAL(UL(J1)))+ABS (AIMAG(UL(J1))) 
J1=J1+N 

IF(RN-ULM) 1,2,2 

RN=ULM 

CONTINUE 

SCL(I)=1./RN 

CONTINUE 

NM1=N-1 

K220 

DO 17 K-1,NM1 

BIG-0. 

DO 11 I=K,N 

IP=IPS(I) 

IPK=IP+K2 
SIZE=(ABS(REAL(UL(IPK)))+ABS(AIMAG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG) 11,11,10 
BIG-SIZE 

IPV-I 

CONTINUE 

IF(IPV-K) 14,15,14 

J-IPS(K) 

IPS (K)-IPS (IPV) 

IPS (IPV)2J 

KPP=IPS(K)+K2 
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16 


17 


PIVOT=UL(KPP) 
KP1=K+1 

DO 16 I=KP1,N 
KP=KPP 
IP=IPS(I)+K2 
EM=-UL(IP)/PIVOT 
UL(IP)=-EM 

DO 16 J=KP1,N 
IP=IP+N 

KP=KP+N 
UL(IP)=UL(IP)+EM*UL(KP) 
CONTINUE 

K2=K2+N 

CONTINUE 

RETURN 

END 


C*****kkkkk*k*FUNCTION BLOG 
C REFERENCE: TECHNICAL REPORT TR-80-1 ;(page 56) 


FUNCTION BLOG(X) 

IFEMET. 71) CO TO 1 

X2=X*X 
BLOG-((.075*X2-.1666667) *X2*1.) *X 
RETURN 

BLOG=ALOG (X+SORT(1.+X*X)) 

RETURN 

END 


QX*kkkkkkkkkkkkkkkkkSGUBROUTINE PLANE 


C 


ເວ 2420) 


REFERENCE: TECHNICAL REPORT TR-80-1 ;(pages 57-64) 


SUBROUTINE PLANE(M1,M2,NP,NT,RH,ZH,XT,AT,THR,R) 


PLANE WAVE EXCITATION VECTOR. FROM MAUTZ AND HARRINGTON. 
MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 


COMPLEX R(1600),U,U1,UA,UB,FA(5),FB(5),F2A,F2B,F1A,F1B 
COMPLEX U2,U3,U4,U5 | 
DIMENSION RH(400),2H(400),XT(4) , AT(4) , R2(4) 

DIMENSION Z2(4),BJ(2000) 

MP-NP-1 

MT-MP-1 

N=MT+MP 

N2=2*N 

CC=COS (THR) 

SS=SIN (THR) 

U=(0.,1.) 

U1=3.141593*U**M1 

M3=M1+1 

M4=M2+3 
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EJ 


23 


20 


19 


Dio 


17 


IF(M1.EQ.0) M3=2 
M5=M1+2 

M6=M2+2 

DO 12 IP=1,MP 

K2=IP 

I=IP+1 
DR=.5*(RH(I)-RH(IP)) 
DZ=.5*(ZH(I)-ZH(IP)) 
D1=SORT (DR*DR+DZ*DZ) 
R1=.25*(RH(I)+RH(IP)) 
IF(R1.EQ.0.) R1-1. 
Z1=.5*(ZH(I)+ZH(IP)) 
DR=.5*DR 

D2=DR/R1 

DO 13 L=1,NT 

R2 (L)=R1+DR*XT(L) 
Z2(L)=Z1+DZ*XT(L) 
CONTINUE 

D3=DR*CC 

D4=-DZ*SS 

D5=D1*CC 

DO 23 M=M3,M4 
FA(M)=0. 

FB(M)=0. 

CONTINUE 

DO 15 L=1,NT 

X=SS*R2 (L) 
IF(X.GT..5E-7) GO TO 19 
DO 20 M=M3,M4 

BJ (M)20. 

CONTINUE 

BJ (2)=1. 

S=1. 

GO TO 18 
M=2.8*X+14.-2./X 
IF(X.LT..5) M=11.8+ALOG10(X) 
IF(M.LT.M4) M=M4 

BJ (M)=0. 

JM=M-1 

BJ (JM)=1. 

DO 16 J=4,M 

J2=JM 

JM=JM-1 

J1=JM-1 
BJ(JM)=J1/X*BJ(J2) -BJ (JM+2) 
CONTINUE 

S=0. 

IF(M.LE.4) GO TO 24 
DO 17 J=4,M,2 

S=S+BJ (J) 

CONTINUE 
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24 
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29 
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26 


2 1 


22 


27 
12 


S-BJ(2)42.*S 
ARG=Z2(L)*CC 
UA=AT(L)/S*CMPLX(COS(ARG),SIN(ARG)) 
UB=XT (L) *UA 

DO 25 M=M3,M4 

FA (M) =BJ (M) *UA+FA (M) 
FB(M) =BJ (M) *UB+FB(M) 
CONTINUE 

CONTINUE 

IF(MI.NE.O) GO TO 26 
FA(1)=FA(3) 
FB(1)=FB(3) 

UA-U1 

DO 27 M-M5,M6 

M7=M-1 

M8=M+1 
F2A=UA*(FA(M8)+FA(M7)) 
F2B=UA*(FB(M8)+FB(M7)) 
UB=U*UA 
FlA=UB*(FA(M8)-FA (M7)) 
F1B=UB*(FB(M8)-FB(M7)) 
U4=D4*UA 
U2=D3*F1A+U4*FA(M) 
U3=D3*F1B+U4*FB(M) 
U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 
R(K2+MT)=-D5*(F2A+D2*F2B) 
R(K5+MT)=D1*(F1A+D2*F1B) 
IF(IP.EQ.1) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.EQ.MP) GO TO 22 
R(K2)=U2+U3 
R(K5)=U4+U5 

K2=K2+N2 

UA=UB 

CONTINUE 

CONTINUE 

RETURN 

END 


C*************SUBROUTINE ZLOAD. 
C REFERENCE: COMPUTATION OF RADIATION AND SCATTERING FOR 


O OO OO 


LOADED BODIES OF REVOLUTION 
HARRINGTON AND MAUTZ 

REPORT AFCRL-70-0046 

AIRFORCE CAMBRIDGE RESEARCH LABS 
CONTRACT NO. F-19628-68-C-0180 
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(309) 2462. 


SUBROUTINE ZLOAD(NP,RH,ZH,Z0,Z) 


COMPUTES IMPEDANCE MATRIX ELEMENTS FOR LOADED BODIES OF REV 
ZO(I) IS THE SURF IMPEDANCE OF THE ITH SEGMENT (NP-1 
SEGMENTS) Z(.) ARE THE IMPEDANCE MATRIX TERMS (TRIDIAGONAL 
FOR T-T SUBMATRIX; DIAGONAL FOR P-P SUBMATRIX). STORED IN 
COL VECTOR. 


10 


41 


42 


GOMEBEXCIT,€2,720(400),/2(2400),X1,X2,X3,V1,V2,X3,FN(400) 
COMPLEX U1,U2,U3,XI,YI 

DIMENSION RH(400),ZH(400),RS(400),D(400),SV(400) 
PI-3.14159 

MT-NP-2 

MP-NP-1 

N=MT+MP 

DO 10 IP=2,NP 

II=IP-1 

DR=RH(IP)-RH(II) 

DZ=ZH(IP)-ZH(II) 

D(II)=SQRT(DR*DR+DZ*DZ) 

SV(II)=DR/D(II) 

RS(II)=.5*(RH(IP)+RH(II)) 

DS=D(II)*SV(II)/2. 

Q1=RS(II)+DS 

Q2=RS(II)-DS 

FN(II)=1. 

IF((ABS(Q2).GT.1.E-6).AND. (ABS(Q1).GT.1.E-6)) 


* FN(II)-ALOG(Q1/Q2) 


CONTINUE 
LO=MT*3-2 

DO 20 I=1,MP 

Cl-PIXZO(I) 

IF(I.EQ.MP) GO TO 80 

KI=2 

IF(I.EQ.1) KI-1 
IF(I.EQ.MT) KI-3 

II=I+1 

C2-PIXZO(II) 

A=SV(I) 
IF(ABS(A).LT.1.E-6) GO TO 41 
X1=C1*FN(I)/2./A 
X22C1*2./A*(1.-RS(I)*FN(I)/D(I)/A) 
X3=-X2*RS(I)/D(I)/A 

GO TO 42 

CONTINUE 
X12C1/2./RS(I)*D(I) 
X2=(0.,0.) 
X3=C1*D(I)/6./RS(I) 
CONTINUE 
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45 


40 
a 


A-SV(II) 

IF(ABS(A).LT.1.E-6) GO TO 45 
Y1=C2*FN(II)/2./A 
Y2=C2*2./A*(1.-RS(II)*FN(II)/D(II)/A) 
Y3=-Y2*RS(II)/D(II)/A 

GO TO 40 

CONTINUE 
Y1=C2/2./RS(II)*D(II) 
Y2=(0.,0.) 
Y3=C2*D(II)/6./RS(II) 
CONTINUE 


C DEFINE TRIDIAGONAL ELEMENTS FOR T-T SUBMATRIX (STORED IN 
C COLS) (Ul- DIAG; U2- LOWER; U3- UPPER) 


XI=X1+X2+X3 

yl=Y1-—Y2+Y3 

IF(KI.EQ.1) XI=C1/SV(I) 
IF(KI.EQ.3) YI=C2/SV(II) 
Ul=XI+YI 

U2=X1-X3 

U3=v] Yo 

L=2+(I-2)*3 

IF(KI.EQ.1) L=0 

Ll=L+1 

L2=L+2 

L3=L+3 
GONEOMĦ(50,60,70)7ku 

Z (L1)=U1 

Ze Ue 

GO TO 80 

Z(L1)-U3 

Z(L2)-U1 

Z (L3)-U2 

GO TO 80 

Z(L1)=U3 

Z (L2)=U1 
Z(LO+I)=2.*C1*D(I)/RS(I) 
CONTINUE 

RETURN 

END 

SUBROUTINE ZTOT(MT,MP,ZL,Z) 


C ADDS THE SURF IMPEDANCE TERMS TO THE TRIDIAGONAL ELEMENTS OF 


C 

a 
90 
60 
70 
80 
20 

€ 

C THE 

6; 


BOR IMPEDANCE MATRIX 2. 


COMPLEX ZL(2400),2(100000) 
N=MT+MP 

MO=MT*3-2 

DO 100 I=1,MP 
LO=MT*N+(I-1)*N+MT 
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IF(I.EQ.MP) GO TO 80 
KI-2 
IF(I.EQ.1) KI-1 
IF(I.EQ.MT) KI-3 
L2=(I-1)*N+I 
L1-L2-1 
L3=L2+1 
M-243*(I-2) 
IF(KI.EQ.1) M-0 
M1=M+1 
M2=M+2 
M3=M+3 
Go fto (5060770) ,K1 
50 Z(D2)=Z2(L2) *ZD(MN1) 
Z(L3)=Z(L3)+ZL(M2) 
GO TO 80 
60 Z(L1)=Z(L1)+ZL(M1) 
Z(L2)=Z(L2)+ZL(M2) 
2 (1L3)S2(L3)+ZL(M3) 
GO TO 80 
OZ (ha \=Z (LL )4+Z25(M1) 
Z(L2)=Z(L2)+ZL(M2) 
80 Z(LO+I)=Z(LO+I)+ZL(MO+I) 
100 CONTINUE 
RETURN 
END 


C*********kSUBROUTINE OGIVE 


CORO OO) DEGREE CO O) 


(9597003 09 


SUBROUTINE  : OGIVE 

DATE a SEPTEMBER 1991 
PROGRAMMER : R.M. FRANCIS 
REVISED : 26 JANUARY 1992 
COMMENTS 


THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION 
(BOR) IN THE FORM OF AN OGIVE. 

DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 

ZH = Z CO-ORDINATE * 2*PI 

RH - RADIUS *2*PI 


SUBROUTINE  OGIVE(NP,ZH,RH,BASE) 


NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 
ZP = ZPRIME, THE POSITION ON Z WHERE THE RADIUS OF 
CURVATURE STARTS 

BASE BASE RADIUS 

RS RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 


INTEGER I,NP 
REAL RH(400),ZH(400) ,ZP, BASE, RS, ZCOORD, RADIUS, AL 
PI=3.1415926 
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C INPUT THE VARIABLES FOR THE OGTVE/ZEAB2E-ANE 


WRITE(6,*) ’ENTER SURFACE CURVATURE (wavelengths) ’ 
READ(5,*) RS 

WRITE(6,*) 'ENTER ZPRIME,WHERE CURVATURE STARTS (wavelengths) ” 
READ(5,*) ZP 

WRITE(6,*) ‘’ENTER BASE RADIUS (wavelengths) ’ 

READ(5,*) BASE 


C PERFORM CALCULATIONS 
AL=SQRT (2. * BASE*RS-BASE**2) 
ZMAX= AL + ZP 
ANG=ASIN(AL/RS) 
L=ANINT ( (RS*ANG+ZP) *5) 
NP=2*L+1 
WRITE (6,*) “NUMBER OF POINTS IS: “,NP 
DZ= ZMAX/FLOAT(NP-1) 
DO 10 I=1,NP 
ZCOORD- (I-1)*DZ 
ZH(I)=2*PI*ZCOORD 
IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 
ELSE 
RADIUS=SQRT(RS**2-(ZCOORD-ZP)**2)+(BASE-RS) 
ENDIF 
RH(I)=2*PI*RADIUS 
10 CONTINUE 
RETURN 
END 


CXXXXXXXXXSUBROUTINE GENEX. 


C SUBROUTINE : GENEX 

C DATE : 14 JANUARY 1992 
C REVISED : 17 February 1992 
C PROGRAMMER : R.M. FRANCIS 


SUBROUTINE GENEX(NM,NP,NT,NPHI,CNRHO,CNPHI,XT,AT,X,A, 
* XR,AR,XP,AP, SCAN, PHI, ARAD, RH, ZH, R) 


SOURCE ON Z AXIS AT Z=0. 

THE EXCITATION VECTOR IS COMPUTED FOR THE GIVEN R, THETA AND 
PHI COMPONENTS SPECIFIED IN FUNCTION SUBROUTINES 
CIRCTHETA, CIRCRHO, CIRCPHI. 


OO €) 


COMPLEX R(1600),CEXP,PSI,ST1,ST2,ST3,ST4,SP,S1,S2,83,584,S5 
COMPLEX CIRCTHETA, CIRCRHO, CIRCPHI 
DIMENSION RH(400) ,ZH(400) ,XT(4) ,AT(4) ,X(100) ,A(100) 
REAL XP(100),AP(100),XR(50) ,AR(50) , SCAN, PHI, ARAD, RHB, ZHB 
INTEGER NM,NP,NT,NPHI,CNRHO,CNPHI,I,N,MP,MT 
PI-3.141592654 
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BK-2.*PI 
MP-NP-1 
MT-MP-1 
N=MT+MP 
C LIMITS ON PHI INTEGRATION 
P1=(2.*PI-0.)/2. 
P2=(2.*PI+0.)/2. 
DO 30 IP=1,MP 
C QUANTITIES FOR THE FIRST SEGMENT (POSITIVE SLOPE) 
I=IP+1 
II-IP 
DR-RH(I)-RH(II) 
DZ-ZH(I)-ZH(II) 
D1=SQRT (DR*DR+DZ*DZ) 
R1=.5*(RH(I)+RH(II)) 
Z1=.5*(ZH(1)+ZH(11)) 
SVP1=DR/D1 
CVP1=DZ/D1 
V1=ATAN2 (DR,DZ+1.E-5) 
C QUANTITIES FOR THE SECOND SEGMENT (NEGATIVE SLOPE) 
C (SKIP IF LAST SEGMENT) 
I=IP+2 
II=IP+1 
DR=RH (I) -RH(II) 
DZ=ZH(I)-ZH(II) 
D2=SORT (DR*DR+DZ*DZ) 
R2=.5*(RH(I)+RH(II)) 
Z2=.5*(ZH(I)+ZH(II)) 
SVP2=DR/D2 
CVP2=DZ/D2 
V2=ATAN2 (DR, DZ+1.E-5) 
C BEGIN PHI INTEGRATION: R TERMS IN S1 AND S2; THETA TERM IN 
C S3 AND S4; PHI TERM IN S5. 


DO 20 J=1,NPHI 
PH=P1*X(J)+P2 
PSI-CEXP (CMPLX (0. ,-MODEXPH)) 


S1=(0.,0.) 
s2=(0.,0.) 
S3=(0.,0.) 
S4=(0.,0.) 
S5=( 110.) 


IF(IP.LE.MT) THEN 
C t-CURRENT INTEGRATION FOR THE POSITIVE SLOPE 
DO 13 L-1,NT 
TP-D1*XT(L)/2. 
RHB=(R1+TP*SVP1)/BK 


91 


ZHB=(Z1+TP*CVP1)/BK 
TH=ATAN2 (RHB, ZHB+1.E-5) 
CC-COS (V1-TH) 
SS=SIN(V1-TH) 

RR=SQRT (RHB**2+ZHB**2) 


S1=S1+AT(L)*(.5+TP/D1)*CC*CIRCRHO(CNPHI, XP,AP, CNRHO, XR, AR, 
k — ARAD,SCAN, PHI,RHB, ZHB) 


S2=S2+AT(L)*(.5+TP/D1)*SS*CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR, 
* ARAD, SCAN, PHI, RHB, ZHB) 
13 CONTINUE 
S1=S1*D1/2. 
S2=S2*D1/2. 
C t-CURRENT INTEGRATION FOR THE NEGATIVE SLOPE 
DO 14 L=1,NT 
TP-D2*XT(L)/2. 
RHB- (R24 TP*SVP2) /BK 
ZHB= (Z2+TP*CVP2) /BK 
TH=ATAN2 (RHB, ZHB+1.E-5) 
SS=SIN(V2-TH) 
CC=COS (V2-TH) 
RR=SQRT (RHB**2+ZHB**2) 


S3=S3+AT(L) *(.5-TP/D2) *CC*CIRCRHO (CNPHI, XP,AP, CNRHO, XR,AR, 
+  ARAD, SCAN, PHI, RHB, ZHB) 


S4=S4+AT(L) *(.5-TP/D2) *SS*CIRCTHETA (CNPHI, XP, AP, CNRHO, XR, AR, 
+  ARAD, SCAN, PHI, RHB, ZHB) 
14 CONTINUE 
S3=S3*D2/2. 
S4=S4*D2/2. 
ENDIF 
C phi-CURRENT INTEGRATION 
DO 15 L=1,NT 
TP-D1*XT(L)/2. 
RHB=(R1+TP*SVP1)/BK 
ZHB=(Z1+TP*CVP1) /BK 
TH=ATAN2 (RHB, ZHB+1.E-5) 
RR=SQRT (RHB**2+ZHB**2) 
S5-S5-AT(L)*CIRCPHI(CNPHI,XP,AP,CNRHO, XR, AR, 
k — ARAD,SCAN, PHI,RHB, ZHB) /R1* (RHB*BK) 
15 CONTINUE 
S5=S5*D1/2. 
SP-SP-A(J)*PSI*S5 
IF(IP.LE.MT) THEN 
ST1=ST1+A(J)*PSI*S1 
ST2=ST2+A(J)*PSI*S2 
ST3=ST3+A(J)*PSI*S3 
ST4=ST4+A(J)*PSI*S4 
ENDIF 


22 


20 CONTINUE 


C 
È 


30 


Ck 


2 ເງ“) 


(9X 9002909 0) 


COMPONENTS ARE STORED IN A COLUMN VECTOR: VT(1,MT), 
VP(MT+1,N) 

IF(IP.LE.MT) THEN 

R(IP)=(ST1+ST2+ST3+ST4)*P1 

ENDIF 

R(IP+MT)=SP*P1 

CONTINUE 

RETURN 

END 


***XXXXFUNCTION CIRCTHETA. 


FUNCTION : CIRCTHETA 

DATE : 1 OCTOBER 1991 

REVISED : 17 February 1992 

PROGRAMMER : R.M. FRANCIS 

COMMENTS : THIS FUNCTION WILL RIGOROUSLY CALCULATE THE 


ELECTRIC FIELD FOR A CIRCULAR APERTURE. THE FIELD IS 
CALCULATED AT THE BOUNDARY DEFINED BY AN OGIVE. THE APERTURE 
IS LOCATED AT Z = 0, AND IS BORE SIGHTED ON THE Z AXIS. THE 
SUBROUTINE OGIVE IS THE SOURCE OF THE GEOMETRIC DATA 
REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 

ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 


COMPLEX FUNCTION CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR,ARAD,SCAN, 


E 


C 
e 


* PHI, RHB, ZHB) 
INTEGER CNPHI, CNRHO,M,N 
REAL XP(100),AP(100),XR(50),AR(50), 
k ARAD, SCAN, PHI, RHB, ZHB 
LOCAL VARIABLES 
REAL Z,RZ, PHIPRIME, PI,S,C, 


k R,RP,X2,Y1,Y2,K 
COMPLEX J. JK Gl, G2 X1,CONI 6C 
* sumx,sumy,sumz,CIRCTHETA 


RR1=(ARAD+0.)/2. 

RR2=(ARAD-0.)/2. 
PI=3.141592654 
K=6.283185307 

J = (0.0,1.0) 

JK= (0.0,6.283185307) 
CON=(0.0,-15.0) 

CC= CON*RR1 


OUTER LOOP: INTEGRATE OVER PHI...-PI<PHI<PI. 
INNER LOOP: INTEGRATE OVER RHO... O0 «RHO«ANT 
RZ=RHB 
Z=ZHB 


S=RZ/SQRT(RZ**2+Z%*%*2) 
C=Z/SORT(RZ**24+Z**2) 
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sumx-(0.0,0.0) 
sumy-(0.0,0.0) 
sumz-(0.0,0.0) 
DO 50 M-1,CNPHI 
PHIPRIME-PI*XP (M) 
DO 60 N=1,CNRHO 
RP=RR1*XR (N) +RR2 


R=SQRT( ((RZ-RP) **2+Z**2)+(4*RZ*RP*SIN ( (PHI-PHIPRIME) /2) **2) ) 


60 
50 


G1-2 (((K*R)**2)-1-(JK*R))/R**3 

G2=(3+(3*JK*R)-(K*R)**2)/R**5 

X1-JK* ((RP*COS(PHIPRIME)*SIN(SCAN))-R) 

Y1=(RZ*COS (PHI))-(RP*COS(PHIPRIME)) 

X2=Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 
sumx=sumx+(CC*AP(M)*AR(N)*(Gi+(X2%*G2))*CEXP(X1)*RP) 
sumy-sumy- (CC*AP(M)*AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz-sumz-t (CC*AP(M)*AR(N) *kY1*Z*G2*CEXP(X1) *RP) 
CONTINUE 
CONTINUE 

CIRCTHETA=(C*COS(PHI)*sumx)+(C*SIN(PHI)*sumy)-(S*sumz) 

RETURN 

END 


C***kkk*kk**FUNCTION CIRCRHO. 


C FUNCTION : CIRCRHO 
C DATE : 1 OCTOBER 1991 
C REVISED : 17 February 1992 
C PROGRAMMER : R.M. FRANCIS 
COMPLEX FUNCTION CIRC" „XP,AP,CNRHO,XR,AR,ARAD, 
* , BHI7RHE, ZHE) 
INTEGER CNPHI,CNRHC 
REAL XP(100),AP(10 am) ,AR(50), 
* ARAD, SCAN, PHI, tuivo, ZHB 
C LOCAL VARIABLES 
REAL 2, RZ, ,PHIPRIME PILI S ,C, 
* R, RE 4275) MUZE 
COMPLEX J 72K, Gi G2) 41, CON ce. 
* sumx, Sumy, sumz, CIRCRHO 
RR1-(ARAD-0.)/2. 
RR2=(ARAD+0.)/2. 
PI=3.141592654 
K=6.283185307 
J = (0.0,1.0) 
JK= (0.0,6.283185307) 
CON=(0.0,-15.0) 
CC-CON*RR1 
C OUTER LOOP: INTEGRATE OVER PHI. ~m =P- PN PIE 
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C INNER LOOP: INTEGRATE OVER RHO... O <RHO<ANT 


RZ-RHB 
Z=ZHB 
S-RZ/SQRT(RZ**2-Z**2) 
C=Z/SORT(RZ**2+Z**2) 
sumx=(0.0,0.0) 
sumy=(0.0,0.0) 
sumz=(0.0,0.0) 
DO 50 M=1,CNPHI 
PHIPRIME=PI*XP (M) 
DO 60 N=1,CNRHO 
RP=RR1*XR(N)+RR2 


R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 


Gl=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**5 
X1=JK*((RP*COS(PHIPRIME)*SIN(SCAN))-R) 
Y1=(RZ*COS(PHI))-(RP*COS(PHIPRIME)) 
x2=Y1**2 
Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 
sumx=sumxX+ (CC*AP (M) *XAR(N) * (G1+(X2*G2)) *CEXP (X1)*RP) 
sumy=sumy+(CC*AP(M)*AR(N)*Y1*Y2*G2*CEXP(X1)*RP) 
sumz=sumZz+(CC*AP(M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 
60 CONTINUE 
50 CONTINUE 
CIRCRHO=(S*COS(PHI)*sumx)+(S*SIN(PHI)*sumy)+(C*sumz) 
RETURN 
END 


C*******FUNCTION CIRCPHI. 


C 


C 
C 
C 


FUNCTION : CIRCPHI 

DATE : 1 OCTOBER 1991 
REVISED : 17 Feb 1992 
PROGRAMMER : R.M. FRANCIS 


COMPLEX FUNCTION CIRCPHI(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 
* SCAN, PHI, RHB, ZHB) 
INTEGER CNPHI, CNRHO,M,N 
REAL XP(100),AP(100),XR(50),AR(50), 
* ARAD,SCAN,PHI,RHB,ZHB 
LOCAL VARIABLES 
REAL Z,RZ,PHIPRIME,PI,S,C, 


* R RP X2,Y1,Y2 RR1, RR2 
COMPLEX J JK,G1,G2%6X1F CON; CO, 
* sumx,sumy,sumz 


PI=3.141592654 
K=6.283185307 
RR1-(ARAD-0.)/2. 
RR2=(ARAD+0.)/2. 

J = (0.0,1.0) 

JK= (0.0,6.283185307) 
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001ໄຂ ( 0 . 0 , -15. 0) 


CC-CON*RR1 
C OUTER LOOP: INTEGRATE OVER PHI... .—-PI<pHI EL 
C INNER LOOP: INTEGRATE OVER RHO... 0 <RHO<ANT 
RZ=RHB 
L-LHB 


S-RZ/SQRT (RZ**24Z**2) 

C=Z/SQRT (RZ**2+Z**2) 

sumx=(0.0,0.0) 

sumy=(0.0,0.0) 

sumz=(0.0,0.0) 

DO 50 M=1,CNPHI 
PHIPRIME-PI*XP (M) 
DO 60 N=1,CNRHO 
RP=RR1*XR (N) +RR2 


R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 


60 
50 


G1=(((K*R)**2)-1-(JK*R))/R**3 

G2-2(3* (G*JK*R) - (K*R) **2) /R**5 

X1=JK*( (RP*COS(PHIPRIME) *SIN(SCAN) )-R) 

Y1=(RZ*COS(PHI))-(RP*COS (PHIPRIME)) 

X2=Y1**2 

Y2-(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 
sumx=sumx+ (CC*AP(M) *AR(N) * (G1+(X2*G2) ) *CEXP(X1) *RP) 
sumy=sumy+ (CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz-sumz-t(CC*AP(M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 
CONTINUE 
CONTINUE 

CIRCPHI=(-SIN(PHI)*sumx)+(COS (PHI) *sumy) 
RETURN 
END 


C***k******SUBROUTINE TESTSPHERE 


C 
C 
E 
C 
C 


SUBROUTINE : TESTSPHERE 

DATE : 21 FEBRUARY 1992 

PROGRAMMER : R. M. FRANCIS 

COMMENTS : GENERATES A SPHERE, WITH PLOTTED POINTS 


AT EQUAL INTERVALS IN THETA. 


SUBROUTINE TESTSPHERE(NP,ZH,RH,BASE,RS,ZP) 
INTEGER NP,I 

REAL ZH(400),RH(400),BASE,RS,ZP,SPRAD 
PI-3.1415926 

BK=2.*PI 

ZP-0.0 

WRITE(6,*)'ENTER AN EVEN NUMBER OF POINTS (NP):’ 
READ(5,*)NP 

WRITE(6,*) “ENTER SPHERE RADIUS’ 

READ(5,*)SPRAD 
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BASE-SPRAD 
RS-SPRAD 

DO 1241 I-1,NP 
ANGLE=PI*FLOAT(I-1)/(2.*FLOAT(NP-1)) 
ZH(I)=BK*SPRAD*COS (ANGLE) 
RH(I)=BK*SPRAD*SIN(ANGLE) 

CONTINUE 

RETURN 

END 

B. SAMPLE DATA FILE OUTLDBOR 


1241 


*** BOR RADIATION PATTERN FOR A CIRCULAR DISC RADIATOR USING 
GENEX *** 


BOR DIAMETER (WAVELENGTHS) =***** 
NUMBER OF GENERATING POINTS (NP)= 70 


SURFACE RADIUS50.00 ZPRIME 0.00 


NT NPHI 
4 60 

INDEX 2(1) RHO (I) SURF IMPED 
1 50.000 0.000 10000.00 0.00 
2 49.987 1.138 10000.00 0.00 
3 49.948 2276 10000.00 0.00 
4 49.883 3.412 10000.00 0.00 
5 49.793 4.547 10000.00 0.00 
6 49.676 5.679 10000.00 0.00 
7 49.534 6.808 10000.00 0.00 
8 49.366 7.934 10000.00 0.00 
9 49.173 9.056 10000.00 0.00 
10 48.954 10.173 10000.00 0.00 
11 48.710 11.285 10000.00 0.00 
12 48.440 12.390 10000.00 0.00 
13 48.146 13 .490 10000.00 0.00 
14 47.826 14.582 10000.00 0.00 
15 47.482 15.667 10000.00 0.00 
16 47.113 16.744 10000.00 0.00 
17 46.720 17.812 10000.00 0.00 
18 46.302 Lo 571 10000.00 0.00 
19 45.861 19 .920 10000.00 0 .00 
20 45.395 20.959 10000.00 0.00 
21 44.906 21.987 10000.00 0.00 
22 44 .394 23 .003 10000.00 0.00 
23 43.859 24.008 10000.00 0.00 
24 43.301 25.000 10000.00 0.00 
25 42 . 2 25.979 10000.00 0.00 
26 42.119 26.945 10000.00 0.00 
27 41.494 27.897 10000.00 0.00 


oy 


40.848 
40.182 
39.494 
38.786 
38.057 
37.309 
36.542 
353755 
34.951 
34.128 
33.287 
32.429 
31.554 
30.663 
29756 
28.834 
274891 
26.945 
25. 979 
25.000 
24.008 
23.003 
21.987 
2027 
9920 
18.871 
1 2.812 
16.744 
15.667 
14.582 
13.490 
128590 
11 285 
102173 
.056 
.934 
.808 
+6079 
„547 
„412 
2716 
71358 
„000 


O > N Q0 > Q O +J 9 


PHI OF RECEIVER (DEG)= 


28.834 
29. 756 
30.663 
31.554 
32.429 
33.297 
34.128 
34.951 
Jio 759 
36.542 


38.057 


39.494 
40.182 
40.848 
41.494 
42.119 
42.721 
43.301 
43.859 
44.394 
44.906 
45.995 
45.861 
46.302 
46.720 
47.113 
47.482 
47.826 
48.146 
48.440 
48.710 
48.954 
49.173 
49.366 
49.534 
49.676 
49795 
49.883 
49.948 
49.987 
50.000 


10000 
10000 
10000 
10000 


10000. 
10000. 
10000. 
10000. 
10000. 


10000 
10000 


10000. 
10000. 
10000. 
10000. 
10000. 
10000. 
10000. 
10000. 
10000. 


10000 


10000. 
10000. 
10000. 
10000. 


10000 
10000 


10000. 
10000. 
10000. 
10000. 


10000 


10000. 
10000. 
10000. 
10000. 
10000. 
10000. 
10000. 
10000: 
10000. 
10000. 
10000. 


0. 


MAXIMUM FIELD VALUE (V/M)= 
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. 00 0.00 
„00 0.00 
„00 0.00 
„00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
„00 0.00 
„00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 000 
00 0.00 
00 0.00 
. 00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
„00 0.00 
„00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
„00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 0.00 
00 


0.16659EF03 


(VOLTS) 
6.9271 
3 .6864 
2.6516 
4.8350 
6.0401 
4.7465 
1.3824 
2.8075 
5.8878 
6.3663 
3.9668 
0.9093 
6.0107 
8.2390 
6.7229 
2.3749 
RA 
9.6622 
10.7713 
6.6302 
3.5922 
7.2792 
13.9956 
12.6783 
14.4138 
1.4466 
5.0101 
18.2410 
22.8067 
17.8045 
18.0824 
9.1689 
25129 
. 3979 
. 4641 
EN 
. 6252 
14 .5891 
15 9270 
45.8795 
78.0162 
106.099 
136.447 
160.822 
161.358 
166.587 
162.565 
153.157 
143.692 


CO=POL 

(DEG) (DB-REL) 
-104 . 51 =2 7. 62 
106: 76 -33. 10 
-50.40 -35.96 
-14.07 -30.74 
=D L =28. 81 
-—2 9D =3 02941 
-4.09 -41 .62 
=174.09 =35. 47 
-174.09 -20.09 
=] 72.58 -28.36 
217/0731 -32.46 
7.83 -45.26 
6.17 -28.85 
4.51 =26. 12 
7.30 - 27.88 
15.98 =36.92 
ALTOS =32.98 
2170.82 -24.73 
=175.50 “2273 
=160. 18 229.00 
-149.07 =33. 33 
2 1:216 2727.19 
303 2.51 
26.06 22. 7 
12.62 =21.26 
93.02 -41.23 
6515 -30 .44 
7161.85 ES S 
20905268 219 27 
2160.35 -19 .42 
-146 .48 EIS 
3.14 25.19 
28.31 =21. 63 
9.45 -14.78 
29.97 = 10.91 
9531525 = 190527 
32. 13 -14.43 
62.48 =2 1. 15 

- 1 70.35351 - 15. 59 
2157.68 -11.20 
2154.23 ສວ... 
7152-60 3492 
2152571 z1:79 
-147.27 7031 
-145.24 -0.28 
-146.15 0.00 
-144.37 zo 
-144.83 “<. “ວ 
-144.21 -1.28 


(VOLTS) 
0.0000 
1.1253 
0.7957 
0.5577 
0.2076 
012 72 
0.5842 
0.7556 
0.9056 
0.9042 
0.6473 
0.2436 
1.0247 
1.2154 
0.9426 
0.4040 
0.5796 
1.3618 
1.5721 
0.6277 
0.8164 
0.9508 
2 1212 
1.1997 
2.0480 
0.5246 
0.7135 
2.3460 
2.5909 
2.0096 
DASS 
1.7552 
0.9059 
4.3773 
5.2557 
3.3511 
212704 
2.0676 
2.8516 
5.2496 
8.1824 

10.9266 

13.7615 

14.3460 

13.9190 

14.2015 

13.4435 

12.5463 

11.3471 


(DEG) 

-92.95 
131 .93 
163 .87 


2505 
2139.66 


5.44 
30.18 
297 
82.84 
.87 
29 


=13 0752 


-96.02 
295309 
. 30 
229 
222 
39 
81.57 
«25 
1257.53 
+15 


-107.94 


271.27 
-86.02 
oo 
2/9 
73 

9308 

80.80 
124.03 


=103.63 
-103.47 
-120. 67 


-84.81 
DD 
37 


X-POL 


(DB-REL) 


zu 
743: 
-46. 
-49. 
-58. 
OS. 
-49. 
-46. 
-45. 
-45. 
-48. 
- GE 
-44. 
TAZE 
-44 . 
2925 


-49 


-41. 
-40. 
-48. 
-46. 
-44. 
dii 
-42 . 
— 36h 
_ 5 0. 
— 4 RE 
= 3 7 
29065 
= 304 
2365 
= 39. 
-45. 
= ele. 
= 0% 
= Je 
-34. 
238: 
mo. 
ES 
Žo. 
-23. 
<21 . 
Zl. 
=z. 
sl. 
=21. 
22. 
23 


65.22 108.369 -140.82 m 
66.52 85.039939. 06 -5.84 
67.83 67. 12 66 E. 82 
69.13 59.199 E45. 8H -8.90 
70.43 47.3499 -145.20 -10.93 
71.74 30.9244 -145.07 -14.63 
73.04 19.4698 =-]43.32 EIS.65 
74.35 12.7001 -149.46 222.96 
75.65 5.20290 - 153 08 EU. 
76.96. 8:/052M71016 5. 2222 
78.26 12.4865 173.28 222750 
79.57 11.0732 91.40 -23.54 
80.87 3.73M PREE =32.94 
82.17 6.6444 103.97 =Ž7.99 
83.48 6.3822 =d 19 23.33 
84.78 3.023393 155.2. -34.82 
86.09 1.8402 75.64 -39.14 
87.39 3.97 0545150. LE 232 2 
88.: /0 7.36099 110705 27:09 
90:00 7.6199 -5.04 226779 


C. SOURCE CODE FOR GENPLT.M 


$ PROGRAM : 


$DATE 


$ PROGRAMMER: 


Z Comments 


A? AV AV 


load bang.m 
load bxpole.m 
load bcpole.m 
load ctest 
for i-1:np 

1f bcpole(i)<=0.0 
bcpole(i)=eps; 


end 
if bxpole(i)<=0.0 


genplt.m 
15 February 1992 


R.M. FRANCIS 


bxpole(i)=eps; 
end 
el(i)210*1og(bcpole(i) 
ep(i)=10*log(bxpole(i) 
zl(i)-10*10g(z(i)); 
lf 21(1 5560 
Z1(1)=-50; 
end 


J; 
) 3 


100 


7.6142 
5.6224 
4.2600 
3.8263 
3.1370 
2.1385 
1.2247 
0.8089 
0.8572 
0.7979 
1.4228 
1.4678 
1.0531 
1.0604 
1.0668 
0.2984 
0.2566 
0.8277 
121926 
1:2270 


65.73 
60.07 
58:69 
67:63 
62.89 
43.54 
3725 
32.28 
=15 °, 
22. 54 
77.84 
-11.88 
= 12 93057 
30252 
-]129 09/0 
104.94 
=5 220 
166.07 
2252 
-108.00 


— 20530 
23443 
31.84 
32-78 
-34.50 
=3 7/83 
-42.67 
-46.28 
-45.77 
=46. 39 
-41.37 
QU TO 
-43.98 
=4 3.92 
-43.87 
-54.94 
=56: 25 
-46.08 
-42.90 
-42.66 


: This program is a plotting routine for the data 
points generated by the main program LDBORMM.F and the 

program TEST.M. Enter parameters iaw the runs of test.m and 
LDBORMM. F. 


if el(i)«-50 


el(i)=-50; 
end 
if ep(i)«-50 
ep(i)—-50; 
end 


end 

susolom(221),plot(deg,z,bang;bepole,'.') 
title('TEST OF LDBORMM.F'),xlabel('degrees') 
Vlabel('magnitude') 

subplot (222) ,plot(deg,zl,bang,el,’.’) 
title(’LOG PLOT (db)’),xlabel(’degrees’) 
ylabel('magnitude') 


text(.1,.4,’analytical (solid), LDBORMM.F(.)’,’sc’) 
text(.1,.3,’PARAMETERS:’,’sc’) 

text(.1,.25,'number of points = 90’,’sc’) 
text(.5,.25,’scan angle = 60 deg.’,’sc’) 
text(.5,.15,’CNRHO = 20’,’sc’) 

text(.5,.2,'CNPHI = 60’,’sc’) 

Eextt(.t,.2, NPHI — 100',;'sec') 

mexe(71,.,15,7NEF = 10” ,7ວເ“) 

text(.l,.1,’/đistance = 50’,’sc’) 
text(.5,.1,'aperture radius 5 57 ,7867) 
text(.1,.05,’complex impedance = (10000.0, 0.0)’,’sc’) 
pause 

aprint 

ແນດ 


load etscat.m 

load epscat.m 

load etf.m 

load epf.m 
subplot(221),plot(bang,etscat) 
Eltile(ETSCAT') 

subplot(222) ,plot(bang,epscat) 
title('EESCAT') 
subplot(223),plot(bang,etf) 
title(ETE/) 
subplot(224),plot(bang,epf) 
title('EPF') 

print 
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APPENDIX B 
TABLE B.1. VARIABLES IN SUBROUTINE GENEXX 


VARIABLE | TYPE DESCRIPTION RANGE 


integer Number of modes in ġ 1< NM < œ 

integer Even number cf points on 2< NP <400 
surface cf the BOR 

integer Even number of points cf 2< c 
integration in T on BOR 

integer Even number of points of 2«NPHI«100 
integration in $ on BOR 


H 


JD D [x a wipe. E, 
Idae | ad QIS 


NPHI 


CNRHO integer Even number of points of 2«CNRHO«50 
integration in p'(source) 

integer Even number of points of 2<CNPHI< 105 
integration in @¢' (source) 

real Abscissa of the function (X) A 

in T,@ ©! and @' 


CRE: 


XT,X,XR, 


AT,A,AR, | real 


Weights of the series 


O<w.<1 
terms(for x,) approximating 
the integral 





SCAN real angular displacement (in -2<4<20 


degrees) of source boresight 
and Z axis 


real orientation in $ of receiver 0«$«180 


ARAD | real dimensionless radius of (paseo 
source antenna (a/À) 





າໆ 
L 
H 





RENE) real 0“ L“ < ae 


3< 1 < 399 


the product of the wave 
number (k) and the radius of 
the BOR, units in radians 





ZH(i) real the product of the wave 0< Zz «o 
number (k) and the distance z | 3< 1 < 399 
on the BOR axis, units in 
radians 


E complex the excitation vector - o« R(i)«e 
1, =2NP-3 4< 1 < 1680 
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APPENDIX C 


A. PROGRAM TO GENERATE DATA FILES FOR GAUSS QUADRATURE 


INTEGRATION 
C PROGRAM : GAUS.F 
C DATE : 12 SEBREMBER 1991 
C PROGRAMMER + R. FRANCIS, D. JENN 
C PROGRAM TO COMPUTE WEIGHTS AND ABSCISSAS FOR GAUSSIAN INT. 


10 


E 


12 


CHARACTER*7 FNAME 
DIMENSION X(500),A(500) 
WRITE(6,*) 'ENTER FILE NAME (gaus---)' 
READ(5,*) FNAME 
WRITE(6,*) ‘ENTER AN EVEN NUMBER (max pts is 500) ! 
READ(5,*) N 
OPEN(2,FILE-FNAME) 
WRITE(2,*) N 
CALL GAUSLEG(-1.,1.,X,A,N) 
DO 10 I=1,N 
WRITE (A)(I) AI) 
STOP 
END 
SUBROUTINE GAUSLEG(X1,X2,X,W,N) 
IMPLICIT REAL*8 (A-H, ,O-2) 
REAL*4 X1,X2,X(N),W(N) 
PARAMETER (EPS=3.D-14) 
M=(N+1)/2. 
XM=.5D0*(X2+X1) 
XL-.5DO*(X2-X1) 
DO 12 I=1,M 
Z=COS (3.1415926535D0* (I-.25D0)/(N+.5D0) ) 
CONTINUE 
P1=1.DO 
P2=0.DO 
DO 11 J21,N 
P3—P2 
P2=P1 
P1=((2.DO*J-1.D0)*Z*P2-(J-1.D0)*P3)/J 
CONTINUE 
PP=N* (Z*P1-P2)/(Z*Z-1.D0) 
Z1=Z 
Z=Z1-P1/PP 
TE(ABS(Z-21).GT.EPS) co TO 1 
X(I)=XM-XL*Z 
X(N+1-I)=XM+XL*Z 
W(I)=2.DO*XL/((1.DO-Z*Z)*PP*PP) 
W(N+1-I)=W(I) 
CONTINUE 
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RETURN 
END 
B. SAMPLE DATA FILE GENERATED BY GAUS.F 
The following is an example of the external sequential 
data file generated by the program GAUS.F. The first row is an 
even integer (N). The data in rows two through N are the 
abscissa and weights for the series approximation of the 


integral calculation. 


Sample data file: gaus20 
20 


< 3912. 1.76140E-02 
205963972 4.06014E-02 
=0. 912284 6.267 206-02 
20:69 8.3276 2E 02 
20: 74633932 1.01950F5E=01 
-0.636054 0.1185 
-0. 510867 0.131689 
1065913 796 0.142096 
=0. 2271786 0.149173 


27:652065E-02 02152 75 


7.65265E-02 0.152755 
0.227786 0.149173 
0.373706 0.142096 
0.510867 0.131689 
0.636054 0.110195 
0.746332 1.0L938E- 01 
0.239117 8.32767E202 
0.912234 6.26720E-02 
07903972 4.06014E-02 
0 9925129 1. JG140FE=0O2 
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APPENDIX D 


A. ARGUMENTS: CIRCTHETA, CIRCRHO AND CIRCPHI 
TABLE D.1. VARIABLE LIST FOR FUNCTION SUBPROGRAMS. 


— = p-e ສະ e 2 


VARIABLE RANGE 
CNPHI sane 
CNRHO same 
XP, XR same 
AP, AR 


real same as RH(i) of same 
Table B.1 

ZHB real same as ZH(i) of same 

m Table B.1 


Eo; E, and E, all 
calculated át the complex 
point on the BOR numbers 
surface defined by 
RHB,ZHB and PHI 













v — 









circtheta 
camcrho 
circphħi 


complex 












B. TEST PROGRAM: CIRCSUB.F 


The following program was used to test for convergence the 
numerical methods used to calculate the electric field 
components for a uniformly illuminated circular aperture. 
Before executing the main program LDBORMM.F, it is desirable 
to determine the minimum number of terms required (CNPHI and 
CNRHO) for the series approximation of the source field to 


converge to the analytical solution of Equation (3.26). The 
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following programs are provided to test for the number of 
terms required for convergence. 

The programs CIRCSUB.F, TEST.M and SCANPLT.M are executed 
in the same directory. To minimize the number of files in the 
main directory, maintain these programs in a subdirectory 
named "circtest". CIRCSUB.F creates the external sequential 
data files ang.m, etheta.m, erho.m and ephi.m. These data 
files are vectors of length NP, containing the angle (0) and 
the magnitude of E,, E, and E,. The matlab program TEST.M 
calculates the analytic solution for E, in accordance Equation 
(3.26) and writes the data to a file CTEST.MAT. The matlab 
program SCANPLT.M loads the data files and plots the results. 


1. Source Code for CIRCSUB.F 


C TESTPROGRAM : CIRCSUDB: E 

C DATE : 1 OCTOBER 1931 

C REVISED : 12 JANUARY 1992 

C PROGRAMMER : R.M. FRANCIS 

C COMMENTS : THIS SUBROUTINE CALCULATES THE ELECTRIC FIELD 
C FOR À CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE 

C BOUNDARY DEFINED BY AN OGIVE. THE APERTURE IS LOCATED AT 2 
C = 0, AND IS BORE SIGHTED ON THE Z AXIS. THE SUBROUTINE 

C OGIVE IS THE SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC 
C TO PERFORM COMPUTATIONS. 

C ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 


C****k*kk** PROGRAM CIRC 

DECLARE VARIABLES INTEGER NP,NPHI,NRHO,I,M,N 

THE CHARACTER STRING PHIPTS/RHOPTS WILL ACCESS THE 
APPROPIATE FILE WITH THE X'S AND COEEFFICIENTS, TO PERFORM 
INTEGRATION BY GAUSSIAN QUADRATURE. EG. GAUS10 YEILDS, 10 
POINTS;GAUS20, 20PTS, ETC. 


ເ) ເງ 6 0 0 


CHARACTER*15 PHIPTS,RHOPTS,GP,GR 
REAL XP(400),AP(400),XR(100),AR(100),S,C,P,PI,SPRAD, 
* SCAN,Z,RZ,PHI,PHIPRIME,DEG,SA,K,ET,ER,EP, 


106 


«E +“ + ss 


Fa 


—n a ša 


— 


* R,RP,X2,Y1,Y2,ARAD,RH(400),ZH(400),angle(400) 
COMPLEX J,JK,G1,G2,X1,CON,CC, 
* sumx,sumy,sum2, 
* ETHETA(400),ERHO(400) , EPHI (400) 
PI-3.141592654 
ISO SO AO) 
JK- (0.0,6.283185307) 
K-6.283185307 
CON=(0.0,-15.0) 
WRITE(6,*) 'ENTER THE RADIUS OF THE ANTENNA, (in wavelengths) ' 
READ(5,*) ARAD 
ANT-ARAD/2. 
WRITE(6,*)'ENTER PHI,(OBSERVATION) IN DEGREES' 
READ(5,*)P 
PHI=P*PI/180. 
WRITE(6,*)'ENTER SCAN ANGLE IN DEGREES' 
READ(5,*)SA 
SCAN--SA*PI/180. 
WRITE(6,*)'ENTER PHIPTS STRING: "gaus### ",(max 400) ! 
READ(5,*)GP 
PHIPTS='gaus/'//GP 
WRITE (6,*) 'ENTER RHOPTS STRING: "gaus### ",(max 100)' 
READ(5,*)GR 
RHOPTS='gaus/'//GR 
C CALL OGIVE(NP,ZH,RH) 
C make loop to enter test sphere. 
write(6,*)'enter the sphere radius' 
read(5,*)SPRAD 
write(6,*)'enter the number of points' 
read(5,*)NP 
do 888 I=1,NP 
angle(I)=PI*(I-1)/(2*(NP-1)) 
RH(I)=SPRAD*SIN(angle(I))*K 
ZH(I)=SPRAD*COS (angle(I))*K 
write(6,*)angle(I),RH(I),ZH(I) 
888 continue 
OPEN(2,FILE-PHIPTS) 
OPEN(3,FILE-RHOPTS) 
OPEN(4,FILE-'etheta.m') 
OPEN(8,FILE-'erho.m') 
OPEN(9,FILE-'ephi.m') 
OPEN(10,FILE-'ang.m') 
READ(2,*)NPHI 
READ(3,*) NRHO 
DO 20 M-1,NPHI 
READ(2,*,END-20) XP(M) , AP(M) 
20 CONTINUE 
DO 30 N-1,NRHO 
READ(3,*,END=30) XR(N) „AR(N) 
30 CONTINUE 
C OUTER LOOP : SOLVES FOR THE SPHERICAL COMPONENTS OF E AT 
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C EACH POINT ON THE SURFACE OF THE OGIVE. 
C MIDDLE LOOP: INTEGRATE OVER PHI...-PI<PHI<PI. 
C INNER LOOP: INTEGRATE OVER RHO... O <RHO<ANT 
DO 40 I=1,NP 
Z=ZH(I)/K 
RZ-RH(I)/K 
S=RZ/SORT(RZ**2+Z**2) 
= Z/SQRT (RZ**2+Z**2) 
DEG-angle(I)*180/PI 
sumx-(0.0,0.0) 
sumy-(0.0,0.0) 
sumz-(0.0,0.0) 
DO 50 M-1,NPHI 
PHIPRIME-ZPI*XP (M) 
DO 60 N-1,NRHO 
RP-(ANT*XR(N) ) x ANT 


R=SQRT((RZ-RP)**2+(Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 

G1-2 (((K*R)**2)-1-(JK*R))/R**3 

G27 (34 (G*JK*R) - (K*R) **2) /R**5 

X1=JK* ( (RP*COS (PHIPRIME) *SIN (SCAN) ) -R) 

Y1=(RZ*COS (PHI) ) -(RP*COS (PHIPRIME) ) 

X2-Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 

CC=CON*ANT 
sumx=Sumx+ (CC*AP(M) *AR(N) * (G1+ (X2*G2) ) *CEXP(X1) *RP) 
sumy=sumy+ (CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz=sumz+(CC*AP(M)*AR(N)*Y1*Z*G2*CEXP(X1)*RP) 

60 CONTINUE 


DIO CONTINUE 
C CALCULATE THE SPHERICAL COMPONENTS AT THE SURFACE. 


ETHETA(I)=((C*COS(PHI)*sumx)+(C*SIN(PHI)*sumy)-(S*sumz)) 
ERHO(I) =((S*COS(PHI)*sumx)+(S*SIN(PHI)*sumy)+(C*sumz) ) 
EPHI(I) =((-SIN(PHI)*sumx)+(COS(PHI)*sumy)) 
ET=CABS(ETHETA(I)) 
ER=CABS(ERHO(I)) 
EP=CABS(EPHI(I)) 
WRITE(4,*)ET 
WRITE(8,*)ER 
WRITE(9,*)EP 
WRITE(10,*)DEG 
40 CONTINUE 
STOP 
END 
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C SUBROUTINE <: OGIVE 

C DATE : 4 SEPTEMBER 1991 

C PROGRAMMER : R.M. FRANCIS 

C REVISED : 13 JANUARY 1992 

C COMMENTS : 

C THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION 
C (BOR) IN THE FORM OF AN OGIVE. 

C DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 

C ZH = Z CO-ORDINATE * 2*PI 

C RH = RADIUS *2*PI 


C * x k k k k k k k k k k k k k k 


SUBROUTINE OGIVE(NP,ZH,RH) 


C NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 
C ZP = ZPRIME, THE POSITION ON Z WHERE THE RADIUS OF CURVATURE 
C STARTS 
C B = BASE RADIUS 
C RS = RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 
INTEGER I,NP 
REAL RH(400),ZH(400),ZP,BASE,RS,ZCOORD,RADIUS,PI 
PI= 3.141592654 
C INPUT THE VARIABLES FOR THE OGIVE,ZP,B,RS,NP 


WRITE(6,*)'ENTER AN EVEN NUMBER OF POINTS, MAX IS 400' 
READ(5,*)NP 
WRITE(6,*) 'ENTER SURFACE CURVATURE (wavelengths)' 
READ(5,*) RS 
WRITE(6,*) 'ENTER ZPRIME,WHERE CURVATURE STARTS (wavelengths)' 
READ(5,*) ZP 
WRITE(6,*) 'ENTER BASE RADIUS (wavelengths)' 
READ(5,*) BASE 
C PERFORM CALCULATIONS 
ZMAX- SQRT((2*BASE*RS)-BASE**2) + ZP 
DZ= ZMAX/(float(NP)-1.) 
DO 10 I=1,NP 
ZCOORD= (float(I)-1.)*DZ 
ZH(I)=2.*PI*ZCOORD 
IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 
ELSE 
RADIUS=SQRT(RS**2-(ZCOORD-ZP)**2)+(BASE-RS) 
ENDIF 
RH(I)=2.*PI*RADIUS 
10 CONTINUE 
RETURN 
END 
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PROGRAM 
DATE : 
REVISED : 


2. Source Code for TEST.M 


TEST. M 

3 OCTOBER 1992 
15 FEBRUARY 1992 
Program TEST.M-- generate analytic soln. 
aperture. 
$saves the data to "mat file", 
npts=input('enter the number of points: ') 
R-input('enter the distance (R):') 


ctest. 


for circular 


r-input('enter the radius (in wavelengths): ') 
scan=input('enter the scan angle(degrees) ') 


ວ 


scan*pi/180; 


eta- 120*pi; 


k-2 


*pi; 


J=sgrt(-1); 
nmax-2*(npts-1); 
cl-j*k*eta*r^2*exp(-j*k*R)/R; 


9 


Lor 


ang (1+1)=pi*i/nmax; 


1=O:npts-1 


deg(i+1)=90*1i/ (npts-1); 
fac=0.5*cos(ang(i+1)); 
if sin(ang(i+1))-sin(sa)== 0.0 


end 


end 


a(i+1)=eps; 


else 


$ calculate the solution from O to 90 degrees.******* 


a(i-1)-22*pi*r*(sin(ang(i-1))-sin(sa)); 


z(irl1)-fac*abs(cl*bessel(1,a(i-1))/a(i-41)); 


np-npts; 


save ctest deg Z np 


3. Data file CTEST.MAT 


0.0000000e+00 
4.7368421e+00 
9.4736842e+00 
1.4210526e+01 
1.8947368e+01 
2.36842116+01 
2.8421053e+01 
3.3157895639 
3.7894737e+01 
4.2631579e+01 
4.7368421e+01 
5.2105263e+01 
5.6842105e+01 
6.1578947e+01 
6.6G315759erT0] 
7.1052632e;01 
7.5789474e+01 


5S.9217626e+03 
5. 405528400 
5.0941450e+01 
4.1944793e+01 
3.1506384e+01 
2.11083666E702 
1.1933350e+01 
4.7301405e+00 
22675065660 
3.2224083e+00 
4.5450714e+00 
4.7331488e+00 
4.2520406e+00 
3.4672834e+00 
2.6244918e+00 
1.8607 87 Ger 60 
T. 230792 T700 


FEO 


20 


8.0526316e+01 7.3451071e-01 
8.5263158e+01 3.4018017e-01 
9.0000000e+01 2.4512433e-16 

4. Source code for SCANPLT.M 


$ PROGRAM : scanplt.m 
$DATE : 15 February 1992 
Comments : This program is a plotting routine for the data 


$ points generated by the test program CIRCSUB.F and the 
$ program TEST.M. Enter parameters iaw the runs of 
$ TEST.M and CIRCSUB.F. 

load ang.m 

load etheta.m 

load ctest 

subplot(211),plot(deg,z,ang,etheta,'.!) 

title('TEST OF SUBROUTINE CIRCTHETA'),xlabel('degrees') 

ylabel('magnitude') 

ນ. ກເ. ດາເ າ າເວັ“ +%1ເແເ611114 (.)','sc') 
pex Ux PARAMETERS:','sc') 
text(.1,.25,'number of points = 90','sc') 
text(.5,.25,'scan angle = 60 deg.','sc') 
Ecc p s CNRHO = 20'2'sc') 
text(. 5,02 WGNPHI =N60" sc!) 
Ee 1-1. 15, distance = 2500','sc!) 
text(.5,.15,'aperture radius = 5','sc') 
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APPENDIX E. MATH NOTES 


A. PLOTS GENERATED BY SUBROUTINE OGIVE 


Plots generated by the subroutine OGIVE do not have equal 
Spacing in the angle 0. As the angle increases the points are 
spaced more closely. Equation 2.1 defines r(z). It follows 


that 





0 = tan“ | Lun | (E.1) 


as Z approaches zero, 0 approaches 90 degrees. The rate at 


which 0 approaches 90 degrees decreases as z decreases since 


2 
do.  Z . (E.2) 
Zz ae 
clearly demonstrates in the limit 
T dor. 
limit = = 0 (E.3) 


zu 


Where the subroutine OGIVE generates the points RH and ZH, the 


density of plotted values increases with angle 0. 
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B. @PLOTS GENERATED BY SUBROUTINE TESTSPHERE 
In the subroutine TESTSPHERE, the variables RH and ZH are 
defined in terms of a sphere radius (SPRAD) and the angle 0 


The following pseudocode applies 


0 (1) E o CAT ii A 0 
RH(i) » SPRAD * sin(0(i)) (E.4) 
ZH(i) » SPRAD * cos (0(1)) 


where A0 is a constant given by 


and NP is the number of points on the BOR. 


A 0 = ree stare (E.5) 


C. NEAR FIELD EQUATIONS FOR A CIRCULAR APERTURE 

The Equations 3.13 to 3.22 were adapted from Equations (6- 
108a) to (6-108c) on page 283 of reference 5. To implement the 
near field integrals for a circular aperture bounded by a BOR 
the cartesian representation of the equations in reference ?? 
were converted to polar coordinates as follows 
= p! cos ó! 


p! sin ġ' (E.6) 
ds’ = dx' dy' = p' dp' dġ' 


MM 
Il 


Ter 


Where the primed notation denotes a source coordinate. 
Straight forward substitution of Equation (E.6) into the 
equations of reference 5 yeild Equations 3.13 to 3.22. 

The cartesian components of the electric field E(x,y,z) 


are transformed to spherical components by 


Es cos(8)cos(4) cos(B8)sin(ġ) -sin(8)] ¡Ex 
E,| = [sin(8)cos (4) sin(0)sin(ġ) cos(0) E, (E.7) 
E} -81ນ (ຜູ່) cos (>) 0 E 


2 


Cylindrical, polar and rectangular 
coordinate systems 





Figure E.1. Polar Geometry for Conversion of E(x,y,z) to 
Spherical Components,E(R,0,4). 


Figure E.1. illustrates the geonetry. 
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